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Abstract 

A  unified  framework  of  continuum  elasticity,  inelasticity,  damage  mechanics,  and  fragmentation  in  crushable  solid  materials  is 
presented.  A  free  energy  function  accounts  for  thermodynamics  of  elastic  deformation  and  damage,  and  thermodynamically  admissible 
kinetic  relations  are  given  for  inelastic  rates  (i.e.,  irreversible  strain  and  damage  evolution).  The  model  is  further  specialized  to  study 
concrete  subjected  to  ballistic  loading.  Numerical  implementation  proceeds  within  a  finite  element  context  in  which  standard  continuum 
elements  represent  the  intact  solid  and  particle  methods  capture  eroded  material.  The  impact  of  a  metallic,  spherical  projectile  upon  a 
planar  concrete  target  and  the  subsequent  motion  of  the  resulting  cloud  of  concrete  debris  are  simulated.  Favorable  quantitative 
comparisons  are  made  between  the  results  of  simulations  and  experiments  regarding  residual  velocity  of  the  penetrator,  mass  of 
destroyed  material,  and  crater  and  hole  sizes  in  the  target.  The  model  qualitatively  predicts  aspects  of  the  fragment  cloud  observed  in 
high-speed  photographs  of  the  impact  experiment,  including  features  of  the  size  and  velocity  distributions  of  the  fragments.  Additionally, 
two  distinct  methods  are  evaluated  for  quantitatively  characterizing  the  mass  and  velocity  distributions  of  the  debris  field,  with  one 
method  based  upon  a  local  energy  balance  and  the  second  based  upon  global  entropy  maximization.  Finally,  the  model  is  used  to  predict 
distributions  of  fragment  masses  produced  during  impact  crushing  of  a  concrete  sphere,  with  modest  quantitative  agreement  observed 
between  results  of  simulation  and  experiment. 

©  2007  Elsevier  Ltd.  All  rights  reserved. 
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1.  Introduction 

An  understanding  of  dynamic  deformation,  damage 
evolution,  and  fragmentation  of  solid  materials  is  needed  in 
order  to  describe  complex  physical  phenomena  occurring, 
for  example,  in  solid  body  collisions  and  ballistic  impacts. 
In  the  context  of  wartime  environments  or  terrorist  attacks, 
injuries  to  soldiers  and  bystanders  due  to  flying  concrete 
debris  as  a  result  of  violently  explosive  destruction  of 
buildings  and  other  urban  structures  have  been  reported 
[!]■  Characterization  of  the  debris  field  would  enable 
enhancement  of  current  protective  strategies,  for  example 
improvements  in  body  armor  and  guidelines  on  safe  stand¬ 
off  distances  from  buildings  undergoing  detonation  or 
pulverization  [2],  The  defense  industry  also  has  a  need  to 
more  fully  understand  the  material  failure  process  so  that 
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strategies  may  be  improved  to  defeat  fortified  concrete 
structures  such  as  reinforced  walls,  bunkers,  and  road¬ 
blocks.  Finally,  crushing  processes  for  rubble  obtained 
from  demolished  urban  structures  (for  recycling  purposes) 
drive  the  development  of  mathematical  models  enabling  an 
increased  understanding  of  inelastic  deformation,  fracture, 
and  fragmentation  of  concrete-based  materials  [3-5]. 

The  aim  of  the  present  study  is  development  of  a  self- 
consistent  theory  accounting  for  dynamic  deformation, 
damage,  and  fragmentation  mechanisms,  specifically  amen¬ 
able  to  brittle,  crushable  solids  such  as  concrete,  mortar, 
and  cinder  block.  This  theory  enables  simulation  and 
analysis  of  urban  structures  undergoing  ballistic  or 
explosive  loading  scenarios.  These  materials  are  referred 
to  here  as  ‘crushable’  since  they  nominally  contain 
significant  initial  porosity.  For  example,  the  microstructure 
of  concrete  consists  of  a  mixture  of  aggregate  stones, 
typically  granite,  quartz,  or  limestone,  embedded  in  a 
cement  matrix.  The  matrix,  alternatively  referred  to  as 
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‘mortar’,  consists  of  sand,  water,  ash,  and  various  binding 
agents.  Porosities  on  the  order  of  10-20%  in  standardized 
concretes  used  for  urban  construction  are  not  uncommon 
[6-8],  and  the  pores  may  be  air-  or  water-filled  depending 
on  the  preparation  method  and  moisture  of  the  ambient 
environment.  This  initial  porosity  induces  a  pressure 
dependency  in  the  effective  compressive  bulk  modulus, 
with  the  stiffness  of  the  material  increasing  as  the  pores  are 
compacted. 

A  voluminous  literature  exists  describing  mechanical 
properties  and  constitutive  modeling  of  concrete.  A  brief 
survey  of  recent  work  deemed  most  relevant  to  the  present 
research  is  given  here.  Hanchak  et  al.  [6]  conducted 
laboratory  triaxial  tests  and  coupled  pressure-shear  experi¬ 
ments,  in  addition  to  ballistic  perforation  measurements. 
Holmquist  et  al.  [7]  developed  a  constitutive  model  (HJC 
model),  with  the  pressure-volume  response  following  data 
from  [6],  and  a  plasticity  and  damage  model  reminiscent, 
yet  not  identical,  to  those  used  previously  for  metallic 
materials  [9],  with  failure  criteria  based  on  cumulative 
strain  measures  [10].  Grote  et  al.  [11]  used  uniaxial  stress 
tests,  split  Hopkinson  bar  tests,  and  plate  impact  tests  to 
measure  the  rate-dependent  compressive  flow  strength  of 
concrete  and  mortar  at  low,  intermediate,  and  high  strain 
rates,  respectively.  Properties  were  subsequently  used  in 
dynamic  finite  element  simulations  of  plate  impact  of  the 
dual-phase  concrete  [12],  in  which  microstructures  were 
resolved  explicitly,  with  an  extended  Drucker-Prager 
plasticity  theory  used  to  capture  pressure-dependent  yield. 
Bazant  et  al.  [13]  developed  a  model  for  concrete  in  which 
plastic  slip  may  occur  on  a  number  of  microplanes, 
somewhat  analogous  to  the  slip  planes  of  crystal  plasticity 
theory  [14],  This  model  has  been  applied  to  address  a 
number  of  features  of  concrete  behavior  arising  in  impact 
events,  including  high  deformation  rates  [15]  and  finite 
strains  [16].  Other  bounding  surface-based  models  addres¬ 
sing  yielding  or  damage  under  a  variety  of  static  and  cyclic 
stress  state  histories  have  been  formulated  [17,18].  The 
variability  in  mechanical  properties  such  as  flow  stress  and 
fracture  toughness  with  micro  structure  constituents  [19], 
processing  conditions  (e.g.,  environment  and  geographic 
location),  and  age  of  the  material  presents  an  inherent 
difficulty  in  precisely  modeling  the  mechanical  behavior  of 
this  class  of  urban  structural  materials. 

The  theoretical  framework  for  the  behavior  of  crushable 
solids  formulated  here  differs  from  many  existing  concrete 
constitutive  models  in  its  usage  of  a  multiplicative  split  of 
the  deformation  gradient  into  elastic  and  inelastic  compo¬ 
nents  [20]  and  its  emphasized  adherence  to  the  laws  of 
continuum  mechanics  and  thermodynamics,  including  the 
entropy  production  inequality  (i.e.,  second  law  of  thermo¬ 
dynamics).  This  is  not  meant  to  imply  that  existing 
engineering  models  developed  elsewhere  for  concrete 
material  behavior  do  not  satisfy  thermodynamic  principles, 
merely  that  such  principles  are  often  not  considered 
explicitly  in  the  process  of  model  development  and 
parameter  selection.  In  the  present  work,  thermodynami¬ 


cally  consistent  properties  and  evolution  equations  for 
porosity  and  damage  are  formulated  following  the  general 
procedures  of  [21,22],  This  approach,  with  multiplicative 
finite  deformation  kinematics  and  simultaneous  adherence 
to  energy  conservation  and  entropy  production,  has  been 
used  frequently  for  metal  plasticity  [23,24],  However,  its 
use  in  urban  structural  materials  such  as  concrete  has  not 
heretofore  been  emphasized  in  the  literature,  though  some 
analogous  thermodynamic  aspects  have  been  incorporated 
for  modeling  geological  materials  [25].  The  present  model 
also  invokes  the  concept  of  an  internal  state  variable 
representing  damage  in  the  material,  related  to  the 
normalized  density  of  micro-cracks  in  the  substance. 
Similar  approaches,  albeit  with  various  different  ways  of 
relating  continuum  damage  variables  to  microscopic 
damage  entities  or  flaws,  have  been  used  for  some  time  in 
continuum  damage  mechanics  theories  [26-30].  Mathema¬ 
tically  consistent  constitutive  models  for  elastic  and 
inelastic  deformation  mechanisms,  including  phase 
changes,  in  moderately  porous  solids  have  also  been 
developed  elsewhere  [31,32]. 

In  the  numerical  implementation  of  the  model  developed 
here,  elastic  strains  are  assumed  to  remain  small  in  order  to 
permit  computational  efficiency  in  an  explicit  integration 
scheme,  as  is  common  in  hydrocodes  used  to  simulate  finite 
plastic  deformation  of  metals  [33].  While  some  consistency 
with  the  theoretical  formulation  is  lost  in  this  approach,  as 
opposed  to  a  more  costly,  yet  more  rigorous  implicit 
formulation  for  finite  elasticity  [34],  possible  limitations  in 
accuracy  due  to  elastic  nonlinearity  are  thought  to  be  more 
apparent  here  than  would  be  the  case  if  small  elastic 
deformations  were  assumed  from  the  outset  in  both  theory 
and  implementation. 

A  variety  of  numerical  approaches  have  been  undertaken 
to  numerically  simulate  ballistic  impact  and  fragmentation. 
Conventional  finite  element  methods  have  been  used  for 
decades  to  model  penetration,  including  concrete  targets 
[7,16].  Often,  elements  are  deleted  or  eroded  when  some 
failure  criteria  or  maximum  cumulative  strain  is  achieved. 
Fragmentation  has  been  addressed  in  this  context  via  post¬ 
processing  calculations  [35,36].  More  recently,  cohesive 
finite  elements  have  been  used  to  simulate  dynamic 
fragmentation  [37,38].  This  technique  is  naturally  more 
realistic  than  element  deletion  for  modeling  discrete  cracks, 
and  is  thought  to  be  particularly  useful  for  simulations  of 
microstructure-level  fracture  along  grain  boundaries  and 
other  internal  interfaces,  for  example,  where  cohesive 
elements  can  be  inserted  along  weak  links  in  the  micro¬ 
structure  [24,39].  The  discrete  element  method,  whereby 
material  elements  are  connected  by  spring-like  entities,  has 
been  used  to  simulate  dynamic  deformation  and  fracture  of 
concrete  [40].  Eulerian  representations  of  material  behavior 
have  also  been  implemented  to  characterize  fragment 
debris  subsequent  to  ballistic  impact  [41,42].  Such  ap¬ 
proaches  offer  advantages  with  regards  to  addressing  fluid¬ 
like  flow  of  material  at  high  pressures.  Particle  methods 
have  also  been  used  to  address  fracture  and  fragmentation. 
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Silling  and  Askari  [43]  implemented  a  theory  in  which  the 
governing  equations  are  cast  in  integral  form.  Smooth 
particle  hydrodynamics  (SPH)  methods  [44,45]  permit 
simulations  of  very  large  distortions  in  ballistic  impact  in 
a  Lagrangian  setting,  with  local  continuum-type  quantities 
such  as  deformation  rates  depending  upon  relative  particle 
motions  and  smoothing  functions.  Concrete  subjected  to 
explosive  loading  has  been  modeled  with  SPH  [45].  The 
approach  followed  here  is  that  of  Johnson  et  al.  [46,47],  in 
which  continuum  finite  elements  are  converted  to  meshless 
particles  when  specified  erosion  criteria  are  met.  To  this 
end,  the  generalized  particle  algorithm  (GPA)  of  the  EPIC 
2003  code  is  used;  GPA  differs  from  SPH  regarding  the 
choice  of  smoothing  functions  [46].  This  method  is 
computationally  efficient  and  appears  natural  for  modeling 
fragmentation,  as  the  particle  velocities  and  trajectories  can 
be  directly  associated  with  those  of  the  fragments  of 
comminuted  material.  However,  information  regarding 
sizes  of  individual  fragments  is  not  readily  available  from 
standard  SPH  or  GPA  methods,  since  the  mass  of  each 
particle  is  simply  the  nodal  mass,  which  in  turn  depends 
upon  the  discretization.  For  example,  a  uniform  grid  would 
produce  a  uniform  size/mass  distribution  of  particles.  This 
issue  is  addressed  here  by  incorporating  additional  physics 
into  the  constitutive  framework  permitting  calculation  of 
non-uniform  fragment  mass  distributions. 

A  number  of  analytical  methods  have  been  developed  to 
quantitatively  characterize  fragmentation.  Grady  [48] 
suggested  that  the  total  energy  of  a  fragmenting  body 
consists  of  expansion  kinetic  energy  and  the  surface  or 
fracture  energy,  the  latter  an  intrinsic  material  property. 
Energy  minimization  under  variations  of  fragment  dimen¬ 
sion  then  yield  the  nominal  fragment  size  for  spherical 
expansion  of  a  fluid  or  brittle  tensile  fracture  of  a  solid.  A 
more  extensive  treatment  was  later  given  in  [49].  Grady’s 
approach  was  extended  by  Glenn  and  Chudnovsky  [50]  to 
include  stored  elastic  energy  and  then  Johnson  and  Cook 
[35]  to  account  for  3D  stress  states,  though  in  these  latter 
approaches  a  direct  energy  balance  was  used,  as  opposed  to 
energy  minimization.  In  a  finite  element  implementation, 
Johnson  and  Cook  [35]  computed  a  cumulative  fragment 
size  based  on  the  strain  and  strain-rate  history  in  the 
material,  though  the  fragmentation  process  was  perhaps 
unrealistically  assumed  to  commence  in  each  element  from 
the  outset  of  local  deformation,  and  the  fragmentation 
energy  consumed  did  not  enter  the  constitutive  model  for 
the  bulk  material  behavior.  Miller  et  al.  [51]  demonstrated 
the  history  dependence  of  fragment  size  using  a  numerical 
approach  with  cohesive  fracture  and  developed  an  analy¬ 
tical  model  of  fragmentation  for  a  body  with  a  time- 
dependent  stress  history.  Analytical  models  have  also  been 
developed  to  characterize  statistical  distributions  of  frag¬ 
ment  sizes.  These  include  methods  based  on  random 
disintegration  of  bodies  in  one  or  more  dimensions  leading 
to  Poisson-type  statistics,  as  well  as  geometry-based 
approaches  partitioning  areas  or  volumes  in  various  ways 
[52,53].  Entropy  maximization  principles,  by  which  the 


most  chaotic  distributions  are  deemed  most  probable,  have 
also  been  used  to  construct  fragment  statistics  [53] 
including  methods  accounting  for  elastic  energy  and 
damage  [54]  and  rotational  inertia  thought  important  for 
granular  microstructures  [55].  Continuum  micromechanics- 
based  models  in  which  crack  sizes  are  related  to  typical 
fragment  sizes  in  dynamically  fracturing  brittle  materials 
have  also  been  developed  [56]. 

Previous  efforts  toward  modeling  ballistic  fragmentation 
with  Lagrangian  or  Eulerian  hydrocodes  have  often 
focused  on  metallic  targets  [35,41,42],  Implementations 
for  modeling  fragmentation  characteristics  of  concrete  or 
geological  materials  have  also  been  reported  [2,36].  In  the 
present  work,  two  alternative  methods  are  considered  for 
computing  fragment  size  and  velocity  statistics,  both 
compatible  with  the  laws  of  thermodynamics  and  momen¬ 
tum  conservation.  The  first  relies  on  a  local  energy  balance 
similar  to  that  of  [35,48-50],  but  newly  applied  to  crushable 
solids  in  a  method  consistent  with  the  description  of  energy 
in  the  bulk  constitutive  model.  Fracture  energy  is  explicitly 
accounted  for  in  the  constitutive  model,  via  evolution  of  an 
internal  variable  representing  damage  in  the  substance,  and 
fuels  the  fragmentation  process  when  damage  reaches  a 
critical  level.  In  this  way,  energy  consumed  in  fragmenta¬ 
tion  is  not  a  static  intrinsic  property,  but  depends  on  the 
damage  progression  in  the  material.  The  distribution  of 
mass  of  the  fragments  is  then  computed  by  application  of 
the  energy  balance  to  converted  particles  in  the  simulation, 
with  the  velocity  distribution  of  the  fragments  associated 
with  that  of  the  particles.  In  the  second  approach,  a 
method  based  on  entropy  maximization  and  classical 
statistical  physics  is  applied,  following  [53,57].  A  joint 
probability  distribution  function  for  fragment  mass  and 
velocity  is  derived  consistent  with  energy  and  momentum 
conservation,  with  the  global  kinetic  energy  and  trajectory 
of  the  fragment  cloud  computed  from  the  mass-weighted 
average  velocity  of  the  ensemble  of  converted  particles. 

The  remainder  of  this  paper  is  organized  as  follows.  The 
continuum  theory  is  presented,  consisting  of  kinematics, 
thermodynamics,  and  kinetic  relations.  Then  fragmenta¬ 
tion  modeling  is  discussed,  including  methods  based  on  a 
local  energy  balance  or  global  statistical  physics.  Features 
that  differentiate  the  present  theory  from  others  in  the 
literature  are  embedded  in  these  descriptions.  Model 
parameters  for  a  particular  mix  of  concrete  are  given, 
and  the  implementation  in  a  finite  element  setting  with 
particle  dynamics  is  then  described.  The  model  is  first  used 
to  simulate  high-speed  perforation  of  a  concrete  target  by  a 
tungsten  sphere.  Numerical  results  are  interpreted  and 
compared  with  experimental  quantities  and  observations 
from  high-speed  photography  [58].  A  second  simulation  is 
also  performed,  whereby  a  concrete  sphere  is  impacted 
against  a  fixed  plate  at  a  moderate  velocity;  resultant 
fragment  masses  are  computed  using  the  local  energetic 
theory  and  compared  with  experimental  data  from  [3], 

The  following  notation  is  invoked.  Cartesian  coordinates 
are  applied  throughout,  with  summation  implied  over 
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repeated  indices.  Vector  and  tensor  quantities  are  repre¬ 
sented  with  boldface  type,  while  scalars  and  individual 
components  of  tensors  are  written  in  italics.  Juxtaposition 
implies  summation  over  two  repeated  adjacent  indices  (e.g., 
(AB)aA  =  The  scalar  product  of  vectors  is  repre¬ 

sented  by  the  symbol  ‘ •  ’  (e.g.,  a*b  —  aaba).  The  colon 
denotes  contraction  over  repeated  pairs  of  indices  (e.g., 
A  :  B  =  tr(ATB)  =  Ac,bBab,  where  ‘tr’  is  the  trace  operation, 
and  C  :  A  =  CabcdAcd).  Superposed  —  1,  T,  and  ‘  •  ’  denote 
inverse,  transpose,  and  material  time  derivative,  respectively. 
Additional  notation  is  defined  as  it  appears  in  the  text. 

2.  Continuum  modeling 

Aspects  of  the  model  framework  for  describing  deform¬ 
ing  and  fragmenting  crushable  solids  are  presented  here. 
The  kinematic  framework  is  general  in  the  sense  that  it  is 
intended  for  applications  involving  any  solid  material 
exhibiting  similar  deformation  mechanisms. 

The  kinematic  description  begins  with  a  multiplicative 
decomposition  of  the  deformation  gradient  F: 

F  =  0x/0X  =  FeFd,  (1) 

where  x  and  X  denote  spatial  and  reference  coordinates  in 
3D  Euclidean  space,  FE  is  the  recoverable  elastic  deforma¬ 
tion,  and  Fd  is  the  irreversible  deformation  associated  with 
defects  such  as  micro-cracks,  voids,  dislocations,  or  shear 
discontinuities  evolving  within  the  material.  The  spatial 
velocity  gradient  then  follows  directly  from  (1)  as 

L  -  8x/8x  -  FEFE~  \  +  FeFDFd~  1  Fe~  \ .  (2) 

le  ld 


where  f  is  the  body  force  vector  per  unit  spatial  volume,  a  is 
the  Cauchy  stress,  and  div  denotes  divergence  in  the  spatial 
frame.  The  spatial  energy  balance  is  written  in  localized 
form  as 


pe  —  a  :  L  —  div  q  (7) 

with  e  the  internal  energy  per  unit  mass  and  q  the  heat  flux 
vector.  The  second  law  of  thermodynamics  is  stated  as 

a  :  L  -  0_lq*Sx0^p(i^  +  0?/),  (8) 


where  0  is  the  temperature,  tj  is  the  entropy  per  unit  mass, 
and  the  Helmholtz  free  energy  is  i Jj  =  e  —  dtj.  The  spatial 
gradient  operator  is  written  as  0X. 

The  Helmholtz  free  energy,  on  a  per  unit  mass  basis, 
is  assumed  to  exhibit  the  following  general  functional 
dependency: 

<A  =  <A(E  E,<P,0,D),  (9) 


where  D  is  a  scalar  internal  state  variable  representing 
cumulative  damage  in  the  material  typically  occurring  in 
conjunction  with  inelastic  deformation.  Stress-strain  and 
temperature-entropy  relations  are  then  deduced  as  [21,22] 


a  —  FEp 


d'A  uET 

0e®f  ' 


>1 


di// 

00' 


(10) 


The  dissipation  inequality  (8)  becomes 


a  :  L° 


(ID 


where  the  Cauchy  pressure  3 p  =  —  trw.  Assuming  isotropic 
heat  conduction  in  the  spatial  frame, 


q  =  -/c0x0, 


(12) 


The  irreversible  volumetric  deformation  associated  with 
pore  collapse  in  crushable  materials  is  described  by 

(p  =  JD~l-  1,  JD  —  det  Fd,  (3) 

where  ip  is  the  volume  reduction  upon  crushing,  a  positive 
quantity  when  the  volume  is  reduced.  Inelastic  volumetric 
expansion  would  be  captured  by  q><  0;  however,  such  an 
effect  is  not  considered  explicitly  in  the  specific  material 
model  that  follows,  since  any  inelastic  tensile  volumetric 
deformations  are  assumed  here  to  remain  infinitesimal  and 
negligible  for  brittle  solids  undergoing  dynamic  fracture. 
The  inelastic  velocity  gradient  from  (2)  can  be  written  as 

Ld  =  FEFDFDIFE1  =  LD  -i<p(l  +<p)_1 1,  (4) 

where  LD  is  the  deviatoric  inelastic  velocity  gradient  mapped 
to  the  spatial  frame  and  1  is  the  identity  map.  The  elastic 
strain  tensor  and  scalar  measure  of  volumetric  elastic  strain 
in  the  intermediate  configuration  are  defined  as 

2Ee  =  FetFe  —  1,  >9E  =  trEE.  (5) 

Standard  local  balances  of  mass,  linear  and  angular 
momentum  apply: 

pQ  =  pj  =  p  det  F,  divn  +  f  =  px,  a  =  nT,  (6) 


with  k  the  scalar  thermal  conductivity,  and  defining  the 
specific  heat  parameter  c  =  de/dd  =  — 002iA/002,  the  energy 
balance  can  be  written  in  terms  of  temperature  rise  0  as 


pcQ  —  o  :  Ld  —  p 
0i A 


0>A  02<A 

dcp  000<p )  ^ 


_0^ 

dD  0007) 


-0- 


000EE 


:EE 


div(K0x0). 


D 


(13) 


The  framework  is  now  specialized  to  brittle  crushable 
solids  such  as  concrete  and  mortar,  with  deformations 
comprising  F°  associated  with  micro-crack  opening  and 
sliding,  as  well  as  pore  collapse  during  compression. 
Specifically,  in  concrete  materials,  deviatoric  deformation 
(represented  here  in  rate  form  by  L  )  consists  of  micro¬ 
cracking,  frictional  crack  sliding,  rubble  formation,  and 
eventual  granular  flow  [12,15,59].  Such  irreversible  defor¬ 
mation  modes  typically  initiate  at  weak  links  in  the 
microstructure  such  as  aggregate-mortar  interfaces,  but 
may  also  occur  alongside  mortar  fissures  and  aggregate 
cracking  at  large  deformations  [4,19].  Let  co  represent 
the  cumulative  local  micro-cracked  area  per  unit  inter¬ 
mediate  volume,  such  that  D  —  m/mc,  where  o>c  is  a 
material  parameter  denoting  the  maximum  sustainable 
crack  density,  subject  to  the  restriction  0^7)^  1.  This  is 
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a  modeling  assumption  unique  to  the  present  theory, 
providing  a  simple  linear  relationship  between  crack 
density,  a  micromechanical  quantity,  and  damage  D ,  a 
macroscopic  continuum  quantity.  When  the  material  is 
undamaged,  the  crack  density  is  assumed  negligible 
(co  =  0),  while  the  crack  density  achieves  its  maximum 
value  (co  =  coc )  when  catastrophic  failure  of  the  material 
element  occurs,  e.g.,  micro-crack  percolation.  Perhaps 
most  ideally,  the  property  coc  could  be  determined  from 
microscopic  investigations,  including  data  obtained  from 
micromechanical  experiments  and/or  computations,  of 
concrete  microstructures  undergoing  dynamic  fracture 
and  fragmentation.  However,  such  microscopic  data  are 
presently  unavailable  for  the  material  of  study  in  the 
simulations  that  follow.  Instead,  coQ  is  chosen  for  concrete 
following  the  experimental  observation  that  typical  frag¬ 
ments  are  of  a  size  commensurate  with  that  of  the 
aggregate,  as  discussed  in  detail  later  in  Section  4.  Such 
an  approach,  with  micromechanically  inspired  constants 
determined  from  macroscopic  data,  is  common  in  practice 
[29].  Though  the  model  presented  here  suffices  for  the 
present  investigation,  more  realistic  relationships  between 
damage  and  microstructure  could  be  envisioned,  at  the 
expense  of  additional  experiments  needed  to  determine  any 
added  parameters. 

The  specific  free  energy  density  is  postulated  on  a  per 
unit  intermediate  configuration  volume  basis  as 

M  =  K(&e,  (p)9l  +  G(  1  -  D) EE:  E® 

+  r(D)  +  Y  (6),  (14) 

where  the  intermediate  mass  density  p  —  pJE ,  K  is  the 
effective  bulk  modulus,  G  is  the  shear  modulus  of  the 
undamaged  material,  E  =  Ee  —  ($e/3)1  is  the  elastic 
strain  deviator,  T  accounts  for  surface  and  concentrated 
elastic  energy  in  the  vicinity  of  micro-cracks,  and  Y 
describes  the  specific  heat  content.  The  bulk  modulus  takes 
the  particular  form 

K  =  KE(cpL  -  cp)/(2cph) 

+  (K\/2  +  K2Se/ 3  +  K3S2E/4)(cp/cpE),  (15) 

where  KE  is  the  elastic  bulk  modulus  of  the  initially  porous 
material,  the  parameter  cpL  denotes  the  maximum  porosity 
reduction  due  to  compressive  pressure  loading,  and  K\,  K2, 
and  K2  determine  the  pressure-volume  relationship  for 
fully  dense  material  at  cp  —  < pL.  The  micro-crack  energy  per 
unit  intermediate  volume  is  written  as 

-r  =  yco  =  (K2ccocD)/(2Ke),  (16) 

where  y  —  K^/2KE  is  the  surface  energy  of  fracture  [49], 
with  Kc  the  effective  fracture  toughness.  The  negative  sign 
denotes  internal  energy  release  (positive  dissipation)  upon 
fracture.  As  will  be  discussed  later,  the  stored  micro-elastic 
energy  released  in  (16)  is  assumed  to  contribute  to  the 
continuum  energy  balance  (13)  only  until  fragmentation 
commences,  following  which  the  energy  dissipated  is 
converted  to  local  kinetic  energy  of  fragment  expansion. 


The  intermediate  second  Piola-Kirchhoff  stress  is 
defined  by 

S  =  p—^-Fr  =  /'EFE_l(rFE_T.  (17) 

f0Ee 

From  (14)  and  (15),  hydrostatic  and  deviatoric  parts  of  S 
are  then  found  as 

p  =  —  trS/3 

=  —  KESE(cpL  —  cp)/cpL  —  (K\Qe  +  K2$E 

+  K2QE)(p/(pL,  (18) 

S  =  S+pl  =  2(1  -Z>)GEE.  (19) 

The  Cauchy  pressure  is  —3 p  =  tr«T  =  tr(/A-1FESFET). 

Deviatoric  plastic  deformation  follows  from  the  flow 
potential  which  is  equated  here  with  the  effective 
deviatoric  stress  S  =  f  (/)ff  :  cr: 

-D  .80 
L  =  X  —  , 

Off 

0  =  [A(l  -  D)  +  B(p/a0f] 

X  [1  +  Cln(s/£o)]cro  (p/(To>-T),  (20) 

where  3a2  =  2L°:  L°  for  stresses  exceeding  the  elastic 

limit,  ff  is  the  deviatoric  stress,  e  =  ^/(|)D  :  D  (with  D  the 

symmetric  deviatoric  part  of  L),  and  A,  B,  C,  N,  ko,  and  cjq 
are  material  parameters.  T  is  equivalent  to  the  ratio  —p/oo, 
with  p  the  tensile  pressure  at  failure.  Note  that  the 
flow  potential  used  in  (20)  is  identical  to  that  of  the  HJC 
model  [7]. 

In  the  present  implementation,  an  isotropic  inelastic 
response  is  assumed  such  that  the  inelastic  spin  may  be 
neglected,  as  indicated  by  LD  =  LDT.  Furthermore,  note 
that  Ld  &  F  Fd~‘  when  elastic  shape  changes  are  small,  an 
assumption  made  later  in  the  numerical  implementation  to 
facilitate  solution  of  large  problems.  The  null  plastic  spin 
assumption  is  standard  for  isotropic  materials  [60].  With 
this  assumption,  the  flow  rule  in  the  first  of  (20),  written  in 
terms  of  the  effective  deviatoric  stress,  is  of  a  form 
consistent  with  the  historic  plasticity  literature  [61]. 

Porosity  and  damage  evolution  are  controlled  via  the 
kinetic  relations 

0  (P^Pc’P^Pi), 

u(p)/v  o  (Pc  <P<Pl)’ 

D  =  kX(\  -  nD(p / gq))  +  T  l(-p/a0),  (21) 

where  a,  k,  and  no  are  positive  constants,  and  the  bracket 
notation  2(x)  =  x  +  \x\.  The  pressure  at  which  inelastic 
crushing  commences  is  denoted  by  pc,  and  the  locking 
pressure  corresponding  to  cpE  is  denoted  as  pL.  Notice 
from  (21)  that  both  cp  and  D  are  always  positive,  i.e., 
irreversible.  Note  that  although  the  particular  forms  of 
Eqs.  (21)  may  be  new,  these  are  differential  equations  of  a 
rate  form  typically  encountered  in  internal  state  variable 
theories  [33]. 
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The  physical  interpretation  and  motivation  for  (21)  is  as 
follows.  At  compressive  pressures  lower  than  pc,  volu¬ 
metric  deformation  is  assumed  elastic  and  pore  collapse 
does  not  occur.  On  the  other  hand,  when  pressures  exceed 
pL,  all  pores  have  been  compressed  and  the  material  is  fully 
dense.  Finally,  at  pressures  between  crushing  (pc)  and 
locking  (pL),  a  linear  relationship  is  implied  between 
changes  in  pressure  and  pore  compaction.  Damage  in  the 
form  of  micro-cracking  evolves  in  conjunction  with 

inelastic  deformation  X  =  \J (|)L°:  L°,  since  inelastic  de¬ 
formation  is  assumed  to  consist  of  matrix-aggregate 
separation  and  micro-crack  sliding.  Increases  in  tensile 
pressure  (p<  0)  further  exacerbate  the  initiation  and 
growth  of  cracks  in  the  microstructure  [59].  It  is  noted 
that  (21)  has  been  developed  for  describing  dynamic  failure 
events,  such  as  spall  and  fragmentation  upon  ballistic 
impact,  and  may  not  reproduce  details  of  rubble  formation 
and  localized  deformation  modes  such  as  cone  and  shear 
failures  occurring  during  static  unconfined  compression  [4], 
for  example. 

The  model  is  next  evaluated  in  terms  of  thermodynamic 
admissibility.  From  (14)  (16), 


P0Z) 


-  E  -  E 

GE  :  E 


Kcmc\ 
2Ke  )' 


(22) 


~Pj^  =  JE~  W  +  *eS|/(2 cpL)  -  (K\/2 

+  K2r%/3  +  K3yE/4)(yE/cpL)l  (23) 

Returning  to  (11),  now  consider  the  strong  form  of  the 
second  law: 


a  :  L 
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Using  (20)  and  (21), 

„  »D  0$  X  „  „ 
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T  \p/c 0)]^0,  (26) 


meaning  that  the  kinetic  relations  for  plastic  flow  (20)  and 
damage  evolution  (21)  are  consistent  with  the  laws  of 
thermodynamics  in  (7),  (8),  and  (24).  Note  also  that  a 
different  flow  potential  could  be  used  in  (20)  should  a 
more  complex  yield  surface  be  necessary,  for  example 
capturing  anisotropy,  creep,  or  cyclic  loading  [15-18].  The 
energy  dissipated  in  (25)  corresponds  physically  to 
irreversible  processes  associated  with  micro-crack  sliding 
and  frictional  granular  flow  mechanisms,  contributing 
mainly  to  temperature  rise  in  (13)  under  near-adiabatic 
conditions.  So  long  as  (25)  is  satisfied,  thermodynamic 
admissibility  of  the  plasticity  component  of  the  model 


framework  is  maintained.  For  evolution  of  porosity,  the 
following  constraint  emerges  in  terms  of  the  pressure, 
considering  (23),  as  cp^O  from  the  first  of  (21): 

P>  (1  +  yLVt/2  +  K2S e/3  +  K3$2E/4)$l  =  A.  (27) 
J  (pL 

Subsequently,  material  parameters  are  selected  such  that 
(27)  is  satisfied  over  the  range  of  applicability  of  the  model. 

Similarities  between  the  present  approach  and  existing 
models  are  clarified  in  the  discussion  that  follows.  Some 
features  of  the  model  developed  here  are  shared  with  HJC 
model  of  [7].  As  will  be  demonstrated  later,  the  EPIC  (2003 
version)  computational  platform  was  used  for  the  large 
scale  simulations,  and  the  built-in  HJC  model  of  that  code 
was  used  as  a  guide  for  developing  some  aspects  of  the 
present  constitutive  model.  To  this  end,  the  authors  of 
[7,62]  are  commended  for  providing  clear  documentation 
and  references  used  in  the  determination  of  material 
constants  of  the  very  same  composition  of  concrete 
(SAC-7)  of  interest  in  the  simulations  that  follow  in  the 
present  study.  Likewise  the  software  developers  are 
acknowledged  for  supplying  an  efficient,  manageable,  and 
portable  code  that  provided  robust  algorithms  for  element- 
to-particle  conversion  [46,47]  essential  for  capturing  frag¬ 
ment  debris  in  the  simulations  to  be  discussed  later. 
Specifically,  here  the  same  plastic  flow  potential  in  (20), 
and  corresponding  plasticity  parameters  (see  later  Table  1) 
are  used  verbatim  from  [7].  However,  the  porosity  equation 
is  written  in  rate  form  (21)  as  opposed  to  a  monotonic 
algebraic  relationship  between  pressure  and  specific  volume 
used  in  [7].  Furthermore,  the  damage  evolution  equation  is 
written  in  rate  form  in  (21),  as  opposed  to  the  incremental 
form  based  on  cumulative  plastic  and  volumetric  strain 
increments  used  in  [7].  Rate  forms  were  used  here  to  satisfy 
the  thermodynamics  of  (13)  and  (24),  which  are  not 
transparently  compatible  with  kinetic  equations  based  on 
cumulative  strain  increments.  It  is  noted  that  the  experi¬ 
mental  data  given  in  [7]  and  references  therein  [6,63,64] 
were  used  to  develop  these  relationships  and  select 
parameters,  as  will  be  discussed  in  more  detail  in  Section 
4.  For  this  reason,  stress-strain  behaviors  predicted  by  the 
present  theory  and  the  HJC  model  are  similar  for  uniaxial 
stress  states,  especially  prior  to  the  accumulation  of 
significant  damage.  However,  it  is  noted  in  [7]  that  the 
constants  for  plasticity  and  damage  evolution  may  not 
have  been  determined  uniquely  in  the  HJC  model  due  to 
lack  of  experimental  data  over  a  complete  range  of 
pressures  and  strain  rates;  the  same  limitations  apply  here. 
For  this  reason,  and  for  the  different  forms  of  the  evolution 
equations  (20),  the  failure  behavior  predicted  by  the 
present  theory  and  that  of  [7]  will  be  different  for 
arbitrarily  general  stress  histories. 

Other  differences  between  the  present  constitutive  model 
and  that  of  [7]  are  now  listed.  The  HJC  model  relies  on  an 
additive  split  of  the  velocity  gradient  into  elastic  and 
inelastic  parts,  as  opposed  to  the  multiplicative  kinematics 
of  (1).  Further,  in  the  HJC  engineering  model,  no  strain 
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Table  1 

Model  parameters  for  SAC-7  (48  MPa)  concrete 


Constant  Value 

Description 

Po 

2440  kg/m3 

Reference  mass  density3 

K 

1.76W/(mK) 

Thermal  conductivity3 

C 

654  J/(kg  K) 

Specific  heat  capacity3 

&E 

15.9GPa 

Initial  elastic  bulk  modulus3 

A'i 

85.0  GPa 

Compressed  first  order  bulk  modulus 

k2 

-151  GPa 

Compressed  second  order  bulk  modulus 

K, 

208  GPa 

Compressed  third  order  bulk  modulus 

k2 

-171  GPa 

Modified  second  order  bulk  modulus3  (pppL) 

G 

14.9  GPa 

Shear  modulus3 

<Pl 

0.10 

Maximum  porosity  reduction3 

Pc 

0.016 

Threshold  pressure  for  irreversible  crushing3 

Pl 

0.080  GPa 

Pressure  at  maximum  porosity  reduction3 

A 

0.79 

HJC  strength  parameter3 

B 

1.60 

HJC  pressure  parameter3 

C 

0.007 

HJC  strain  rate  parameter3 

N 

0.61 

HJC  pressure  exponent3 

£o 

1.0/s 

Reference  strain  rate3 

<70 

0.048  GPa 

Static  compressive  yield  strength3 

f 

0.083 

Maximum  allowable  tensile  pressure3 

Kc 

0.831  MPa  m1/2 

Fracture  toughness 

Wc 

1 .7(10s)  m  1 

Maximum  crack  density  (coat  D  =  1) 

a 

0.00614 

Porosity  evolution  parameter 

k 

300 

Damage  evolution  parameter 

nD 

4.00 

Damage  evolution  parameter 

Dj 

0.50 

Threshold  damage  for  fragmentation  initiation 

“Denotes  parameter  obtained  or  derived  from  Ref.  [7], 


energy  function  is  given  from  which  the  elastic  moduli  are 
derived.  Instead,  piecewise  linear  and  nonlinear  relations 
are  given  for  the  bulk  modulus,  fit  to  the  data  of  [6], 
Modifications  to  the  shear  modulus  upon  damage  accu¬ 
mulation  or  porosity  compaction  are  not  discussed  in  the 
original  paper  [7].  On  the  other  hand,  here  a  free  energy 
function  is  postulated  explicitly  in  (14),  from  which  the 
elasticity  relations  are  derived  naturally  in  (10),  (15),  and 
(17)  (19).  Here,  the  Green  elastic  strain  (5)  is  used  in  the 
free  energy  (9)  as  an  independent  constitutive  variable, 
consistent  with  other  accepted  theories  of  finite  deforma¬ 
tion  thermo-mechanics  [23,24,34]. 

The  shear  modulus  is  reduced  linearly  with  increasing 
damage  D  in  (14)  and  (19),  inspired  from  simple  damage 
mechanics-based  arguments  relating  effective  moduli  with 
micro-crack  densities  [26-28].  In  the  HJC  theory  [7],  no 
connection  between  D  and  micro-cracks  is  made.  However 
the  assumption  made  here,  relating  D  to  micro-crack 
density,  is  not  necessarily  inconsistent  with  the  yield  surface 
(20)  formulated  in  [7].  Furthermore,  the  prescription  of 
linearly  decreasing  deviatoric  inelastic  strength  in  (20)  with 
increasing  D  is  consistent  with  the  analogous  reduction  in 
shear  modulus  with  increasing  D  in  the  present  theory. 

Also,  in  (15)  the  effective  bulk  modulus  is  interpolated 
linearly  between  elastic  (KE)  value  at  null  compression 
and  the  cubic  Hugoniot  curve  at  full  compaction.  In  the 
present  framework,  the  bulk  modulus  is  not  reduced  upon 
accumulation  of  damage  Z>,  but  such  an  effect  could  be 


incorporated  to  reflect  damage-induced  losses  in  volu¬ 
metric  elastic  stiffness  under  tensile  stress  states.  It  is  noted 
that  similar  pressure-volume  responses  are  predicted 
between  the  present  theory  and  the  HJC  model  [7]  for  a 
particular  composition  of  concrete,  since  both  are  fitted  to 
the  same  triaxial  compression  data  [6]. 

A  fundamental  difference  of  the  present  theory  with 
many  others  [7,12,15-18]  is  the  present  model’s  explicit 
tracking  of  internal  energy  changes  due  to  evolution  of 
damage  D  in  the  constitutive  theory.  The  fracture  energy  T 
in  (14)  and  (16)  accounts  for  the  local  elastic  energy 
released  by  breaking  bonds  during  micro-crack  initiation 
and  extension.  This  energy  is  distinct  from  that  dissipated 
during  plastic  deformation  in  (25),  e.g.,  crack  sliding  and 
granular  flow  contributing  to  frictional  heating.  However, 
both  dissipative  mechanisms  tend  to  occur  simultaneously, 
as  damage  growth  is  driven  by  plastic  flow  via  (21),  except 
for  cases  wherein  the  compressive  pressure  provides 
enough  confinement  to  suppress  increases  in  D  in  the 
constitutive  model.  As  will  be  discussed  later  in  Section  3, 
this  fracture  energy  provides  a  counterbalance  to  the 
expansion  kinetic  energy  that  drives  the  fragmentation 
process  (under  the  same  deformation  conditions,  the 
greater  the  fracture  energy,  the  fewer  fragments  produced 
per  unit  mass).  Such  an  effect,  whereby  the  energy 
dissipated  due  to  damage  evolution  in  the  bulk  constitutive 
model  is  explicitly  transferred  to  fuel  the  formation  of 
fragment  debris,  has  not  been  emphasized  elsewhere.  In 
particular,  the  HJC  model  [7]  does  not  at  present 
incorporate  any  means  for  computing  fragment  sizes, 
though  some  post-processing  capabilities  exist  in  the  EPIC 
code  for  metals  [35].  Also,  no  connection  between  D  and 
crack  energy  is  made  in  the  HJC  model  [7]. 

Tangible  benefits  arise  from  formulating  a  finite  defor¬ 
mation  constitutive  theory  that  explicitly  accounts  for 
restrictions  imposed  by  the  laws  of  thermodynamics,  as 
opposed  to  a  purely  mechanical  model  that  may  appear 
more  simple  and  tractable  for  general  engineering  practice. 
When  selecting  parameters  and  constitutive  equations, 
non-physical  behavior  can  be  avoided  if  thermodynamic 
restrictions  on  dissipation  are  addressed  [65].  Irreversible 
processes  such  as  damage  accumulation,  internal  friction, 
and  inelastic  deformation  should  cause  energy  dissipation, 
as  opposed  to  energy  storage,  in  the  material  [16].  Incorrect 
estimates  of  energy  dissipation  could  lead  to  inaccurate 
temperature  predictions  in  high  rate  simulations  of 
material  behavior,  when  conditions  may  be  nearly  adia¬ 
batic.  Effects  of  temperature  on  the  mechanical  response 
will  be  greatest  if  the  elastic  constants,  pressure,  and/or 
flow  stress  depend  explicitly  on  temperature.  They  do  not 
here,  but  such  temperature  effects  could  be  incorporated 
later  into  the  present  theory  when  data  become  available, 
and  such  effects  are  thought  to  be  important  for  crushable 
materials  subjected  to  shock  loading  [32]. 

The  value  of  the  present  general  theory  may  be  more 
fully  realized  as  it  is  extended  to  other  materials.  For 
example,  finite  deformation  elasticity  may  be  even  more 
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useful  for  describing  those  materials,  such  as  crystalline 
ceramics,  that  undergo  little  plasticity  or  inelastic  volume 
change  under  shock  loading.  The  connection  of  damage 
variable  D  to  micro-crack  density  provides  an  opportunity 
for  computational  and  experimental  micromechanics 
[28,56]  to  provide  a  physical  basis  for  damage  evolution 
equations  and  parameters. 

While  the  theory  appears  complex,  it  is  emphasized  that 
the  inequalities  entering  the  thermodynamic  analysis  in 
(24)-(27)  need  not  be  consulted  once  the  kinetic  equations 
and  corresponding  parameters  have  been  chosen  for  a 
particular  material.  It  is  also  noted  that  the  number  of 
parameters  in  the  present  theory  is  not  excessive,  with  only  a 
few  more  needed  here  than  in  [7]  for  example,  specifically  in 
order  to  address  fracture  energy  and  fragmentation. 
Furthermore,  in  the  large  scale  computations  to  be  discussed 
later,  much  of  the  execution  time  is  apparently  consumed  by 
particle  neighbor  search  algorithms,  as  opposed  to  the 
constitutive  update  of  the  material  response. 

3.  Fragmentation  modeling 

Two  methods  for  describing  fragment  size  and  velocity 
statistics  are  developed  here.  The  first  method,  termed  the 
‘energetic  approach’,  relies  on  a  local  energy  balance  to 
compute  the  fragment  size  and  number  of  fragment(s) 
associated  with  each  local  volume  element  (e.g.,  a  finite 
element,  computational  cell,  or  particle  in  a  numerical 
scheme),  and  assigns  to  all  fragment(s)  the  local  velocity  of 
that  volume  element.  The  second  method,  termed  the 
‘statistical  physics  approach’,  is  derived  upon  maximiza¬ 
tion  of  a  global  statistical  entropy  function  subject  to 
constraints  regarding  conservation  of  mass,  energy,  and 
momentum. 

3.1.  Energetic  approach 

Fragmentation  of  the  local  volume  element  is  assumed  to 
occur  over  a  (small)  time  period,  beginning  when  the 
damage  reaches  a  threshold  value,  D  =  Dj.  Over  this 
period,  the  energy  released  per  unit  volume  due  to  internal 
micro-cracking  is  estimated,  from  (16),  as 

rD 

peD=  /  [K2co)C/(2Kn)\  dD 
Jdt 

=  K2ccoc(D  -  Dt)/(2Ke).  (28) 

During  the  fragmentation  process,  the  material  is  assumed 
to  retain  all  energy  apart  from  cd-  which  may  be  stored  via 
(14)  or  dissipated  as  heat.  Prior  to  fragmentation,  F 
contributes  to  the  dissipation  and  possible  temperature  rise 
via  (13)  and  (26);  then,  upon  D  —  Dj ,  this  energy 
contributes  to  the  relative  kinetic  energy  of  fragments. 
Instantaneous  energy  transfer  from  fracture  to  fragmenta¬ 
tion  follows  from 

MeD  =  Mp~lt\D>DT  =  d^fw!  +  ms)/ dr,  (29) 


where  M  is  the  total  mass  of  the  volume  element,  and  u\ 
and  ms  denote  absolute  energies  due  to  the  relative  linear 
and  spin  momenta  per  fragment.  Note  that  these  energies 
are  distinct  from  the  kinetic  energy  of  translation  and 
rotation  of  the  center  of  mass  of  the  element,  which  are 
assumed  to  be  conserved  during  the  fragmentation  process. 
Summation  in  (29)  is  implied  over  all  fragments  comprising 
the  particular  volume  element  under  consideration. 

Assuming  cube-shaped  fragments  with  edge  length  h, 
energies  can  be  estimated  as 

u\—b5ph2/ 16,  us  —  b5  pf /12,  (30) 

where  the  effective  strain  rate  in  the  fragment  is 
e  =  VdTd,  with  D  the  symmetric  part  of  L,  and  (j)  is  the 
rate  of  rotation  of  the  fragment  about  its  moment  of 
inertia.  The  velocity  gradient  L  is  assumed  to  be  distributed 
uniformly  over  the  fragment  population  comprising  the 
continuum  volume  element,  equivalent  to  the  global  value 
assigned  to  that  element  [35].  Combining  (29)  and  (30), 

(X^3)^D  =  ytY)ib5pk2/l6  +  h5P^/lT>-  (31) 

Assuming  that  the  mass  density  is  equal  and  constant 
among  N  fragments,  replacing  b  with  a  mean  effective 
length  b  —  (pN)~x  M,  and  integrating  with  respect  to  time, 

(31)  becomes 

KeP(s)2 

where  the  rotation  cj)  has  been  omitted  in  (32),  as  the 
translational  energy  is  assumed  here  to  far  exceed  the 
rotational  energy  of  the  fragments,  and  since  rotation  rates 
of  fragments  are  difficult  to  quantify  from  experimental 
observations  [58]  modeled  later.  Equation  (32)  yields  the 
mean  fragment  dimension  b.  Mass  conservation  then 
provides  the  local  number  of  fragments  N=M/(pb  ). 
Note  that  (32)  is  a  time-integrated  generalization  of  the 
energy  rate  balance  (31);  as  a  result,  the  mean  fragment 
dimension  is  inversely  proportional  to  the  strain  rate  in 

(32) ,  and  not  the  second  time  derivative  of  strain. 
Equation  (32)  may  be  written  in  the  alternative  form 

b  =  A[K2c/(KEps2)f,  (33) 

where  A  —  2.83y/a>c(D  —  Dj)  and  N  —  /.  This  can  be 
compared  with  the  brittle  fragmentation  model  of  Grady 
[48],  in  which  A  =  2.71  and  N  =  I,  and  the  ductile  failure 
model  of  Johnson  and  Cook  [35],  which  exhibits  the  form 
of  (33)  with  the  substitutions  A  =  4,  N  =  and 
K2c  =  Ke6e.  All  models  predict  a  decrease  in  mean 
fragment  dimension  with  increasing  strain  rate,  in  agree¬ 
ment  with  general  observations  on  dynamically  fragment¬ 
ing  solids  [49].  Decreasing  fragment  sizes  with  increasing 
strain  rate  have  been  predicted  elsewhere  specifically  for 
concrete  and  mortar  [66]. 

Note  that  the  derivation  of  (32)  follows  from  the 
assumption  that  the  expansion  kinetic  energy  of  the 
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fragments  balances  the  fracture  energy  of  the  material; 
similar  assumptions  have  been  made  elsewhere  with  minor 
variations,  most  notably  [35,48-51].  The  theory  leads  to  the 
logical  conclusion  that  the  more  rapidly  a  mass  of  material 
is  straining  or  deforming,  the  more  fragments  it  will 
produce.  The  use  of  effective  strain  rates  generalizes 
derivations  based  on  spherical  expansion  or  ID  failure 
[48-50].  The  angular  rotation  rate  of  the  fragments  is 
neglected  in  the  applications  that  follow  in  the  present 
work,  though  such  rotations  could  logically  derive  from  the 
spin  of  the  continuum  material.  For  example,  analogously 
to  the  definition  of  e,  one  may  assume  that  (j)  =  |e  :  W|, 
where  W  is  the  skew  part  of  L,  and  e  is  the  rank  3 
permutation  tensor. 

The  new  contribution  here  is  introduction  of  parameters 
Dj  and  a>c  into  the  fragment  size  expression  (32).  The 
former  introduces  history  effects  in  a  simple  manner:  only 
the  strain  rate  e  that  occurs  during  the  fragmentation 
process,  when  D>Dj,  is  used  in  the  computation  of  the 
fragment  size,  as  is  discussed  in  more  detail  in  Section  4. 
This  prevents  fragments  from  forming  during  processes  in 
which  the  material  undergoes  no  damage,  for  example 
during  purely  elastic  deformation.  The  use  of  a>c  connects 
the  fragmentation  energy  consumed  during  dynamic 
expansion  with  the  fracture  energy  resulting  from  damage 
evolution  in  the  bulk  constitutive  model. 

A  common  approach  involves  only  post-processing 
calculations  to  determine  fragment  sizes,  with  the  average 
strain  rate  over  the  entire  history  of  the  problem  considered 
[35].  In  such  an  approach,  the  fragmentation  energy 
does  not  influence  the  constitutive  response  of  the  bulk 
material.  On  the  other  hand,  here  the  fragmentation 
event  influences  the  bulk  response  and  vice  versa  through 
the  exchange  of  energy,  from  damage  dissipation  to 
fragment  expansion,  as  discussed  following  Eq.  (28) 
and  illustrated  in  Fig.  1.  In  general,  theories  such  as 
[35,48-51]  are  based  on  an  energy  balance  or  energy 
minimization  and  are  consistent  with  thermodynamic 
principles.  However,  implementation  of  such  models 
in  a  purely  post-processing  capacity,  irrespective  of 
how  the  fragmentation  energy  is  supplied  or  extracted 
from  the  bulk  constitutive  model,  may  not,  in  principle, 
properly  account  for  the  energy  consumed  during  frag¬ 


mentation.  It  is  noted,  however,  that  such  energy  ex¬ 
changed  during  fragmentation  may  be  negligible  compared 
to  that  induced  by  plastic  dissipation,  and  that  material 
elements  undergoing  significant  damage  accumulation 
tend  to  support  little  mechanical  strength,  so  such  a  non¬ 
physical  assumption  may  present  few  difficulties  in 
practical  computations. 

3.2.  Statistical  physics  approach 

Here,  the  entire  fragmenting  body  is  considered  at  once, 
and  it  remains  to  be  determined  how  the  masses  and 
velocities  are  distributed  among  fragments  comprising  this 
body.  The  issue  is  resolved  upon  consideration  of  entropy 
maximization  constrained  by  mass,  momentum,  and 
energy  conservation.  In  what  follows,  fragment  mass  and 
velocity  distributions  are  derived  individually,  then  com¬ 
bined  to  form  a  joint  probability  function  capturing  both 
fragment  sizes  and  speeds. 

The  mass  distribution  follows  from  entropy  maximiza¬ 
tion  subjected  to  the  constraints  that  the  total  mass  and 
total  number  of  fragments  are  known,  following  a 
procedure  outlined  by  Grady  and  Kipp  [53].  These 
constraints  are  written  as 

M  =  p  ]T  b3  =  pNb\  (34) 

where  it  is  assumed  that  the  density  is  the  same  among  all 
fragments  and  the  velocity  does  not  affect  the  mass 
distribution,  and  M  and  h  are  known  from  mass  conserva¬ 
tion  and  a  global  application  of  the  preceding  energetic 
analysis,  respectively.  A  measure  of  the  statistical  entropy 
associated  with  the  mass  distribution  is  Sm  =  A: a  In  P , 
where  k b  is  Boltzmann’s  constant  and  P  is  the  number  of 
possible  fragment  arrangements.  The  value  of  Sm  is 
maximized  by  maximizing  the  quantity  [67] 

J 

InP  =  Ain N  —  ^  n,  Inn,-,  (35) 

i=0 

where  n,  is  the  number  of  fragments  of  mass 
m,  ^  m  <  mi  +  dmh  with  5ml  describing  the  range  of  masses 
admitted  in  each  bin  i,  and  with  j  the  total  number  of  bins. 
In  (35)  it  is  implicitly  assumed  that  Hi  are  large  numbers. 


Undamaged 
D  =  0 


Damaged 
0  <D<Dt 

crack  energy  and  friction 
— *  temperature  rise 


Fragmenting 
D<  D  <  1 

crack  energy  release 

— » fragment  expansion  KE 


Fig.  1.  Damage  evolution  and  fragmentation  process. 
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Constraints  (34)  are  rewritten  simply  as 

j  i 

N  ='^2  nu  M  —  m‘- 


(36) 


i= 0 


i= 0 


Introducing  Lagrange  multipliers  a  and  /?,  an  equivalent 
function  /(«,•)  is  defined  as 

f  =  \nP  +  oi(N-Yjn)+p(M-Yjm(),  (37) 

which  is  maximized  via  the  solution  of 


=  -(1  +  In rii)  -  a  -  /ira,  =  0. 
C«; 


(38) 


Making  the  notation  change  a  +  1  — >■  a,  (38)  yields 

itj  —  exp(— a  —  pmj).  (39) 

Rewriting  (36)  as  integrals  of  continuous  functions  [53], 


N=  n  dm  =  /  exp(— a  —  /fin)  dm 

Jo  Jo 

=  exp(-a)/0, 


(40) 


p  OO  /*  oo 

M  =  /  mndm=  /  m  exp  (—a  —  pm)  dm 

Jo  Jo 

=  exp(— a)//)2.  (41) 

Finally,  (39)  becomes 

«(/«,)  =  n,  =  (N2/M)exp(—Nmj/M),  (42) 

with  the  cumulative  probability  distribution  of  fragments 
larger  than  m  given  by 


h(m)/N  =  No  exp(-Nom), 


(43) 


where  NqM  =  N. 

The  fragment  velocity  distribution  follows  from  similar 
arguments.  A  measure  of  the  statistical  entropy  associated 
with  the  velocity  distribution  is  Sy  =  k b  In  W,  attaining  its 
greatest  value  upon  maximization  of  the  function 


In  W  —  N  In  N  —  n>  ln 


(44) 


/=o 


where  n,-  is  the  number  of  fragments  with  kinetic  energy 
e,  <  e  <  c,  +  Seh  with  Se,  spanning  the  range  of  energies 
admitted  in  each  bin  i.  The  constraints  on  the  distribution 
are  written  as 


n = ^2  n>’  E = X!  e‘ ’ 


(45) 


1=0 


;=0 


where  E  is  the  total  kinetic  energy  of  the  fragment  cloud 
that  will  be  determined  later.  Following  an  analogous 
procedure  to  that  in  (37)  (38), 

nj  =  exp(— a  -  /(<?,),  (46) 

where  the  Lagrange  multipliers  i  and  are  determined  as 
follows.  Making  the  analogy  with  rigid  molecules  [67],  let 

[}=  \/(kRT)  =  (3N)/(2E),  (47) 


where  T  is  a  thermodynamic  temperature.  Then 

«,■  =  exp[-a  —  (3Nei)/(2E)\.  (48) 

Assume  that  in  ballistic  scenarios  the  fragment  velocity  is 
unidirectional,  coaxial  with  the  velocity  of  the  center  of 
mass  of  the  fragment  cloud,  such  that  2e,-  =  Then  the 
probability  distribution  is 

h{vi)  —  A  exp[(—  3Nmitf)/(4E)\,  (49) 

where  A  —  exp(— a)  is  determined  by  normalization: 

/>oo  poo 

1=  /  h{v)dv  =  A  /  exp[(— 3Nmv2)/(4E)\  dvt 

Jo  Jo 

=  A \J  (jiE)/ (3mN),  (50) 


giving 


h(v)  —  \J  (3m  N) /( nE )  exp[(—  3  Nmv2)  /  (4  E)\ . 


(51) 


The  form  (51)  differs  from  the  velocity  distribution 
obtained  by  Grady  and  Winfree  [57],  who  assumed  that  the 
fragments  may  scatter  randomly  in  three  dimensions, 
yielding  3D  Maxwell-Boltzmann-type  velocity  statistics. 
In  the  ballistic  concrete  perforation  experiment  that  is 
described  and  simulated  later,  the  flying  debris  tend  to 
follow  a  roughly  uni-directional  path,  with  fragments 
enclosed  in  a  cone  angle  of  less  than  45°,  as  opposed  to  a 
path  of  fully  radial  expansion. 

The  joint  probability  distribution  of  mass  and  velocity  is 
derived  by  combining  (43)  and  (51): 


p(m,  v)  =  h(m)n(v) 


/3mA5 

nEM 2 


exp 


N 


3  N 
4 El 


(52) 


The  following  relations  then  arise  for  total  probability, 
total  fragment  mass,  mean  velocity,  and  total  linear 
momentum: 


poo  poo 

/  /  p(m,  v)dmdv  =  1, 

Jo  Jo 

poo  poo 

/  /  /}(/?/,  v)m  dm  dv 


(53) 


Jo  Jo 


j  mn(m)yj  ii(v)dv  ]  dm  =  M, 


(54) 


N . 


p(m,v)vdmdv 


o  Jo 


l  r 0°  /  r°°  \  4 e 

=  nJ„  ”(0)1’dljdm  =  Vw' 


(55) 


p(m,v)mv  dmdv 


Jo  Jo 


>  /  r°°  \ 

mh(m)  (  /  n(v)v  dv  \  dm  —  JEM/3 .  (56) 

Identifying  (56)  with  the  linear  momentum  of  the  fragment 
distribution,  conservation  of  linear  momentum  demands 
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that  the  energy  distributed  to  feed  the  velocity  probability 
distribution  is 

E  =  7>pb^Nv2.  (57) 

The  mean  velocity  in  (55)  then  becomes  2v,  i.e.,  twice  the 
velocity  of  the  center  of  mass  of  the  fragmenting  body.  The 
expansion  energy  of  (28)  does  not  contribute  to  the  velocity 
distribution  because  this  energy  is  consumed  by  the  strain 
rate  of  the  fragments  as  given  by  Eq.  (30),  as  opposed  to 
their  linear  velocities. 

The  method  presented  above,  based  on  global  entropy 
maximization,  is  offered  as  an  alternative  to  the  local 
energy  balance  method  derived  in  Section  3.1.  The  second 
method  provides  a  basis  for  comparison  of  numerical 
results  in  the  absence  of  quantitative  data  on  fragment 
velocities,  which  are  scarce  for  materials  of  present  interest. 
The  differences  between  the  two  approaches  are  as  follows. 
For  the  approach  in  Section  3.1,  a  local  energy  balance 
leads  to  a  local  fragment  dimension.  In  the  numerical 
implementation  to  be  discussed  later  in  Section  4,  each 
particle  is  then  assigned  its  own  number  of  fragments,  with 
the  size  of  these  based  on  the  local  strain  rate  history  and 
damage  evolution  in  that  particle  (or  corresponding 
element).  Consideration  of  all  particles  then  provides 
statistical  distributions  of  fragment  mass  and  velocity, 
with  velocities  of  fragments  obtained  directly  from  parent 
particles.  In  contrast,  in  the  statistical  method  of  the 
present  section,  distributions  of  fragment  sizes  and 
velocities  are  found  using  (52),  irrespective  of  the  local 
particles’  velocities  and  the  masses  associated  with  each; 
only  the  total  fragment  mass  and  kinetic  energy  are  needed 
(these  are  supplied  by  the  computation).  In  this  way, 
distributional  information  is  provided  by  the  universal 
entropy  maximization  procedure  in  lieu  of  the  local  physics 
from  the  constitutive  update  and  particle  algorithms.  This 
tends  to  make  the  statistics-based  method  less  sensitive  to 
the  meshes  used  in  the  computations,  as  will  be  demon¬ 
strated  in  Section  5.1,  though  the  sensitivity  of  the  former 
energetic  approach  to  mesh  density  is  not  excessive  here. 
Both  methods  are  motivated  by  and  address  principles  of 
mass,  energy,  and  momentum  conservation.  The  new 
contribution  in  the  present  section  is  derivation  of  the 
joint  mass-velocity  distribution  (52)  in  the  context  of  ID 
fragment  trajectories,  and  derivation  of  Eq.  (57),  the  energy 
fueling  the  fragment  velocities,  neither  of  which  was  given 
explicitly  in  [57]. 

It  is  noted  that  geometry-based  fragmentation  theories, 
such  as  discussed  in  [52,53],  are  often  not  developed  via 
direct  thermodynamic  considerations.  However,  this  does 
not  mean  such  models  contradict  thermodynamic  princi¬ 
ples,  and  in  fact,  many  have  been  analyzed  and  justified 
using  entropy  methods  from  statistical  mechanics  [54—57]. 
So  long  as  the  energy  consumed  during  damage  and 
fragmentation  is  tracked  correctly  in  the  constitutive 
model,  it  is  thought  here  that  such  approaches  are 
reasonable  from  a  thermodynamic  perspective.  In  the 


problems  considered  later  in  Section  5,  calculations 
accessing  the  constitutive  model  are  needed  to  supply 
global  data — specifically  the  collective  mass,  number,  and 
linear  momentum — entering  the  statistical  theory  to 
produce  the  distributions  of  mass  and  velocity,  which  in 
turn  satisfy  global  conservation  laws  via  (53)  (57).  Though 
unlikely,  if  such  data  are  known  a  priori  from  experimental 
or  analytical  means,  then  the  constitutive  model  and 
numerical  simulations  will  not  be  needed  for  predicting 
fragmentation  statistics  with  a  single  equation  such  as  (52) 
or  those  in  [52,53]. 

4.  Concrete  material  modeling 

The  model  was  applied  here  to  study  a  particular 
concrete  material.  Properties,  parameter  selection,  and 
numerical  implementation  are  described  briefly  in  what 
follows. 

4.1.  Material  properties  and  parameters 

The  preceding  theory  was  applied  to  describe  a  concrete 
of  unconfined  compressive  strength  48  MPa  (7  ksi),  as 
studied  previously  by  Hanchak  et  al.  [6]  and  Holmquist  et 
al.  [7],  Parameters  are  listed  in  Table  1. 

Constants  marked  with  a  superscript  were  obtained  from 
[7];  most  of  these  enter  the  yield  function  (20).  Note  that 
our  framework  is  general  in  the  sense  that  an  alternative — 
and  possibly  more  robust  for  certain  complex  load 
histories — plastic  potential  for  deviatoric  flow  could  be 
substituted  for  and  the  HJC  plasticity  parameters 
entering  Eq.  (20)  while  still  maintaining  thermodynamic 
admissibility  of  the  kinetic  relations.  Also,  identical 
constants  were  often  found  to  adequately  address  behavior 
for  which  elastic  strains  were  small,  such  that  differences  in 
pressures  p  and  p  in  intermediate  and  spatial  frames, 
respectively,  were  negligible. 

The  fully  crushed  material  is  assumed  to  behave  similarly 
to  the  aggregate  at  high  pressures,  following  [7].  For  this 
concrete,  the  aggregate  is  assumed  to  consist  of  fine-  and 
coarse-grained  granite  stones  [6].  Constants  K \ ,  Ki,  and 
K},  describing  the  volumetric  elastic  response  of  the  fully 
crushed  material,  were  determined  by  fitting  the  pressur¬ 
e-volume  response  of  (18)  to  the  shock  Hugoniot  data 
for  granite  [64],  assuming  purely  volumetric  elastic 
deformation  of  the  form  FE  =  (1  —  ;)\.  The  fit,  shown  in 
Fig.  2,  is  deemed  accurate  for  pressures  up  to  20GPa. 
Porosity  evolution  parameter  a  was  found  from  linear 
interpolation  between  crushing  and  locking  pressures,  i.e., 
a  =  (Pl^o/IPl  ~  Pc)-  The  pressure-volume  response  under 
triaxial  loading  at  pressures  under  20GPa  is  shown  in 
Fig.  3  and  compared  with  experimental  data  [6].  Here,  //  = 
P/Po  —  1  is  the  total  volumetric  strain.  The  complete 
pressure-volume  response  and  porosity  change  under 
triaxial  compression  are  given  in  Fig.  4;  notice  the  abrupt 
increase  in  effective  bulk  modulus  upon  attainment  of  the 
locking  pressure,  at  which  cp  =  cpL.  Thermodynamic 
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Fig.  2.  Pressure  p  vs.  elastic  volumetric  deformation  x,  fully  crushed 
material  (experimental  data  from  [64]). 


Fig.  5.  Thermodynamic  dissipation  ratio  associated  with  porosity 
evolution. 


Fig.  3.  Pressure  p  vs.  volumetric  deformation  p,  small  compression  range 
(experimental  data  from  [3]). 


Fig.  4.  Porosity  compaction  <p  and  pressure  p  vs.  volumetric  deformation 
p,  large  compression  range  (model  only). 

admissibility  of  porosity  evolution,  from  inequality  (27), 
was  verified  for  all  tp> 0,  as  shown  in  Fig.  5. 

Fracture  toughness  Kq  was  obtained  from  [67]  and  was 
assumed  constant.  However,  it  is  noted  that  variability  in 
fracture  strength  is  expected  depending  upon  loading  rate, 
stress  state,  specimen  size,  and  composition  [68].  Following 
monotonic  and  cyclic  slow-rate  test  data  and  analyses  given 
in  [6,7,63],  failure  occurs  ( D  =  1)  when  the  cumulative 


plastic  strain  1  =  0.0033  under  null  pressure  conditions 
p/a o  =  0,  giving  k  =  300.  Consulting  the  above  references, 
prescribing  D  =  1  when  ).  =  0.01  under  average  pressure 
p/a0  =  i  yields  nD  =  4. 

The  computed  stress-strain  responses  for  the  concrete 
material  deformed  in  uniaxial  compression,  tension,  and 
shear  are  shown  in  Fig.  6.  Results  were  obtained  from 
deformation  of  a  single  dynamic  finite  element  under 
conditions  where  imposed  components  of  the  deformation 
gradient  F  and  nominal  strain  e  were  related  as  follows:  for 
compression  Fp  =l-e,  tension  F\ \  =  1  +  g,  and  shear 
F 12  =  £.  The  axial  stress  is  plotted  for  tension/compression 
and  the  shear  stress  for  shear  deformation.  Note  that  the 
material  hardens  slightly  with  strain  rate,  particularly  in 
compression,  in  agreement  with  recent  experiments  [11]. 
The  rate  effect  is  less  noticeable  in  pure  shear,  as  the  curves 
for  rates  e  of  1/s  and  10/s  are  nearly  superposed.  Softening 
and  complete  stress  relaxation  occur  under  tensile  and 
shear  loading  due  to  damage  accumulation.  Slight  oscilla¬ 
tions  occur  in  the  tension  and  compression  stress-strain 
curves  due  to  lateral  stress  relief  waves  from  the  Poisson 
effect  and  the  dynamic  integration  scheme. 

Maximum  crack  density  coc  and  threshold  damage 
parameter  Dj  were  chosen  based  on  observations  from 
ballistic  experiments  on  the  material  of  interest  [58].  In 
experiments,  the  typical  fragment  size  was  observed  to  be 
on  the  order  of  the  minimum  dimension  of  the  coarse 
aggregate  of  the  concrete  microstructure,  here  9.5  mm. 
For  cubic  fragments,  this  implies  an  edge  length  b  of 
9.5/V3  mm,  as  the  aggregate  stones  were  sized  in  practice 
via  passage  through  a  sieve  with  holes  of  diameter  9.5  mm. 
A  typical  strain  rate  observed  over  the  duration  of  the 
fragmentation  event  was  s  =  2(  10)4/s,  based  on  examina¬ 
tion  of  results  from  simulations  reported  later  in 
Section  5.1.  Invoking  (32),  and  assuming  that  fragmenta¬ 
tion  begins  when  the  material  has  lost  half  of  its  strength, 
i.e.,  Dj  =  0.5,  then  produces  the  value  of  n>c  listed  in 
Table  1.  More  quantitative  experiments  on  damage  evolu¬ 
tion  and  fragmentation  are  needed  to  support  independent 
choices  of  co c  and  D  \  :  however,  for  the  present  scenarios  of 
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Fig.  6.  Computed  dynamic  stress-strain  (a  —  e)  response. 


interest,  the  choice  of  parameters  insignificantly  affects  the 
model  predictions  apart  from  the  mean  fragment  dimen¬ 
sion.  Others  have  suggested  a  correlation  between  grain  or 
aggregate  size  and  typical  or  minimum  fragment  dimension 
in  granular  geological  materials  [69], 

4.2.  Numerical  implementation 

The  material  model  was  inserted  into  the  2003  version  of 
the  EPIC  Lagrangian  finite  element  code  [47,62].  The 
equations  of  motion  and  constitutive  response  are  inte¬ 
grated  explicitly,  and  GPA  and  contact  algorithms  are 
available  for  addressing  multi-body  interactions. 

The  stress  rate  is  derived  by  differentiating  (17)  and 
assuming  small  elastic  stretch,  giving 

6  =  2G(  1  -  T>)DE  +  WEo  -  oWE  -  Da/(  1  -  D),  (58) 

where  D  is  the  symmetric,  deviatoric  part  of  LE  and  WE  is 
the  skew  part  of  LE.  The  assumption  of  small  elastic  strains 
is  standard  in  metal  plasticity  literature  [33],  has  been 
widely  used  for  concrete  modeling  [12],  and  enables 
efficient  integration  of  the  deviatoric  stress  state  via  a 
radial  return  algorithm.  The  assumption  is  thought  to  be 
justified  here  for  concrete  deforming  in  shear  or  tension  in 
that  yielding  and  failure  of  the  material  should  occur  prior 
to  the  attainment  of  large  deviatoric  stresses  (and 
correspondingly  large  deviatoric  elastic  strains),  from 
evolving  inelastic  deformation  (20)  and  damage  D  in  (21). 

Under  realistic  cases  where  elastic  volume  changes  are 
large,  e.g.,  triaxial  compression,  the  hydrostatic  pressure  is 
of  primary  interest.  At  large  pressures,  the  hydrostatic 
response  is  integrated  in  terms  of  the  modified  total 
volumetric  strain  fi,  following  [7]: 

P  =  (li~  <PL)/(1  +  <Pl)-  (59) 

At  pressures  beneath  the  locking  pressure,  elastic  strains 
are  small,  and  a  relationship  equivalent  to  (18)  is  used.  On 
the  other  hand,  a  direct  cubic  fit  to  the  Hugoniot  data  of 
[64]  is  invoked,  via  p  =  K\fi  + K2p2  +  ,  at  high 

pressures  and  large  elastic  compressive  strains,  after  all 
porosity  is  compacted.  In  this  way,  possible  errors  due  to 


linearization  of  elastic  strain  propagated  from  (58)  are 
avoided  in  the  high  pressure-volume  response.  Note  from 
Table  1  that  AT  differs  from  AT  by  20GPa  due  to  the 
choices  of  volumetric  strain  measures  used  at  low  (5)  and 
high  (59)  pressures. 

The  energy  balance  (13)  is  exercised,  along  with 
thermodynamic  expressions  (23)  and  (26),  to  update  the 
temperature  of  the  material.  To  this  end,  the  elastic  strain 
quantities  of  (5),  and  hence  the  elastic  deformation 
gradient  FE,  are  needed.  The  latter  is  computed  via  an 
exponential  update  standard  in  the  context  of  computa¬ 
tional  crystal  plasticity  [24,70].  The  following  approxima¬ 
tions  are  used  here: 

Ff+A,  =  exp(LEA?)Ff , 

exp(LEAf)  =  1  +  S'^~  LEA?  +  - — C°S^(Af)2(LE)2,  (60) 

where  At  is  the  time  increment  of  integration,  LE 
is  assumed  constant  over  the  time  increment,  and 

Q  —  \J (1  /2)Le:  LeAa  In  the  present  application  of  (60), 
Le  is  computed  by  subtracting  the  plastic  deformation  rate 
from  the  total  velocity  gradient,  thereby  including  all  rigid 
body  rotations  in  the  elastic  part  of  the  deformation.  One 
can  show  that  the  form  of  (60)2  is  exact  when  LE  is  skew, 
and  otherwise  introduces  numerical  error  of  order  (A/)3 
that  is  negligible  in  the  context  of  dynamic  calculations  of 
the  sort  here,  for  which  At  is  very  small.  Also,  since  the 
elastic  stretch  generally  remains  small  in  the  large  scale 
computations  that  follow  in  Section  5,  in  these  computa¬ 
tions  Le  consists  predominantly  of  rotation  and  is  nearly 
skew. 

The  following  methodology  for  addressing  failure  of  the 
material  is  used.  When  an  integration  point  achieves  a 
critical  value  of  damage,  i.e.  when  D  =  1,  failure  occurs. 
Failed  material  supports  no  deviator  stresses  or  tensile 
pressures.  In  practice,  when  D  >  0.95  in  the  numerical 
integration,  D  is  subsequently  set  to  unity  and  the  material 
fails  in  the  next  integration  cycle,  in  order  to  avoid 
difficulties  with  (58)  as  D  approaches  unity.  Additionally, 
finite  elements  are  converted  into  particle  nodes  when  a 
scalar  measure  of  effective  inelastic  or  volumetric  strain, 
termed  the  erosion  strain,  is  attained  [47,62],  in  order  to 
alleviate  numerical  inaccuracies  associated  with  highly 
distorted  elements.  In  the  calculations  that  follow,  the 
erosion  strain  is  chosen  as  0.5,  following  recommendations 
in  [62].  For  the  present  material  model,  failure  may  or  may 
not  necessarily  precede  erosion,  but  eroded  elements  are 
almost  always  failed,  as  D  in  (20)  typically  attains  a  value 
of  unity  prior  to  accumulation  of  sufficient  strain  needed  to 
trigger  particle  conversion. 

Two  methods  have  been  developed  for  addressing 
fragment  mass  and  velocity  distributions.  For  the  approach 
of  Section  3.1,  a  typical  fragment  dimension  is  computed 
for  each  converted  particle  using  Eq.  (32),  where  the 
strain  rate  s  during  the  fragmentation  event  is  time- 
averaged  only  over  explicit  integration  cycles  for  which 


282 


J.D.  Clayton  /  International  Journal  of  Impact  Engineering  35  (2008)  269-289 


D  >Dj.  The  number  of  fragments  associated  with  a  given 
particle  is  then  found  by  dividing  the  particle’s  mass  by  the 
mass  of  its  associated  fragment.  Each  fragment  is  assigned 
the  velocity  of  its  parent  particle.  For  the  statistical 
approach  of  Section  3.2,  the  total  mass  and  mass-weighted 
velocity  of  the  particle  cloud  (see  later  Eq.  (62))  are 
extracted  from  the  simulation  output.  The  total  number  of 
fragments  follows  from  assumption  of  a  nominal  fragment 
dimension  corresponding  to  the  aggregate  size  as  discussed 
above.  Then  upon  application  of  (57),  all  quantities  in 
distribution  (52)  are  known,  and  the  velocity  and  mass 
probability  distributions  are  computed  in  a  post-processing 
step. 

The  presentation  here  has  emphasized  the  general 
theory,  with  more  restrictive  assumptions  that  narrow  the 
focus  of  the  model  to  a  particular  concrete  material 
introduced  only  as  they  become  necessary.  This  enables 
avenues  for  extending  the  framework  to  related  materials 
in  the  future,  and  also  permits  description  of  the  limitations 
of  the  specific  material  model  and  its  numerical  imple¬ 
mentation.  Note  that  while  the  nonlinear  elastic  theory  was 
implemented  numerically  for  hydrostatic  deformations  in 
order  to  select  bulk  elastic  constants  (Figs.  2  and  3)  and 
judge  the  thermodynamic  admissibility  of  porosity  evolu¬ 
tion  (Fig.  5),  some  simplifications  have  been  enacted  to 
enable  solution  of  large  scale  problems  encompassing 


Table  2 

Steps  used  in  constitutive  update  for  large  scale  computations 


Step 

Equation(s) 

Read  in  constants,  internal  variables,  and 
rate  variables 

- 

Compute  effective  shear  modulus 

G  -»  G(1  -D) 

Compute  effective  flow  stress 

(20) 

Compute  deviatoric  Cauchy  stress 

a  =  a  —  (tro/3)l 

Update  stress  using  radial  return 
algorithm 

(58) 

Update  inelastic  deformation  rates 

(2),  with 

2Ld  =  FDFd-‘  +  FD“rFDr 

Update  plastic  work 

(25) 

If  ip  =  q> L,  compute  effective  bulk 
modulus  and  pressure 

(15),  (18) 

If  q>  =  tp L,  compute  effective  volume 

(59),  with 

change  and  pressure 

p  =  K\ji  +  K2p.2  +  K3p? 

Update  elastic  deformation  gradient 
Compute  elastic  volume  change  and 
elastic  strain  deviator 

(60) 

(5),  EE  =  Ee  -  (,9e/3)1 

Update  porosity  and  damage 

(21) 

Add  dissipation  from  rates  of  porosity 
and  damage 

(23),  (26) 

Update  internal  energy 

(13) 

If  D>Dt.  update  cumulative  strain  during 
fragmentation 

e/+A/  =  £/  +  At^/D  :  D 

If  D>Dt ,  update  rate  over  fragmentation 
time  span,  t— tp 

E=  £,+A ,/(t+At-  tF) 

If  D>Dt.  compute  fragment  size  and 
number  of  fragments 

(32),  N  =  M/(ptJ‘'\ 

If  D  =  1 ,  zero  deviator  stresses  and  tensile 
pressure 

<t  ->  0,  p  e  (p< 0)  ->  0 

hundreds  of  thousands  of  elements  and  particles.  Specifi¬ 
cally,  these  include  Eqs.  (58),  (59),  and  the  failure 
algorithm  discussed  above.  Table  2  lists  the  general  scheme 
and  corresponding  equations  used  for  the  constitutive 
update  in  the  large  scale  calculations  presented  in  Section 
5.  These  steps  in  the  computation  are  bypassed,  for  the 
most  part,  once  local  failure  of  the  material  has  occurred. 
In  Table  2,  t  =  ty  denotes  the  time  instant  in  the  solution  at 
which  fragmentation  commences,  when  the  damage  D 
attains  or  surpasses  a  threshold  value  of  Dj.  It  is  also  noted 
that  all  numerical  implementations  rely  on  some  approx¬ 
imation  of  the  governing  differential  equations,  e.g., 
discretization  in  space  and  time,  and  the  present  computa¬ 
tions,  in  the  context  of  coupled  finite  elements  and  meshless 
particles,  are  of  course  subject  to  such  universal  limitations. 

5.  Numerical  simulations 

Two  problems  were  used  to  demonstrate  the  capabilities 
of  the  concrete  material  model.  The  first  involved  perfora¬ 
tion  of  a  concrete  target  by  a  spherical  metallic  projectile; 
the  second  involved  fragmentation  of  a  spherical  concrete 
projectile  upon  impact  with  a  metallic  plate. 

5.1.  Concrete  perforation 

A  specific  initial-boundary  value  problem  was  solved  to 
mimic  a  ballistic  perforation  experiment  conducted  at  the 
US  Army  Research  Faboratory  (ARE)  [58].  In  the 
experiment,  a  small  tungsten  sphere  was  fired  at  concrete 
wall.  The  sphere,  of  diameter  7.94  mm,  was  formed  from  an 
alloy  of  mass  density  18  690  kg/m3.  The  concrete  target  was 
SAC-7  composition  [6],  with  no  reinforcing  bars,  25.4  mm 
thick.  High-speed  photography  was  used  to  record  the 
striking  velocity  of  the  projectile,  measure  the  residual 
velocity  of  the  penetrator,  and  collect  images  of  the 
fragment  debris  cloud.  Post-mortem  measurements  of  hole 
and  crater  dimensions  were  obtained  via  visual  inspection. 

The  initial  problem  geometry  is  shown  in  Fig.  7. 
The  concrete  target’s  dimensions  were  25.4  mmx 
102  mm  x  102  mm,  with  half  modeled  explicitly  as  a  result 
of  symmetry.  Two  meshes  were  used  to  evaluate  effects  of 
mesh  density:  a  coarse  grid  with  102  144  composite 
tetrahedral  elements  [62]  and  a  fine  grid  with  244  512 
elements.  The  mesh  density  in  the  target  was  reduced 
gradually  with  increasing  distance  from  the  location  of 
impact.  The  impact  velocity  of  the  sphere  was  1120  m/s. 
Frictionless  contact  between  projectile  and  target  was 
assumed.  Adiabatic  conditions  were  invoked,  with  a 
uniform  initial  temperature  of  294  K. 

The  constitutive  model  used  for  the  concrete,  as 
developed  in  Sections  2-4,  is  the  focus  of  this  work  and 
was  used  to  model  the  target.  On  the  other  hand,  a  simple 
yet  efficient  and  widely  used  constitutive  model  was  chosen 
for  modeling  the  behavior  of  the  projectile.  Specifically,  the 
alloy  of  the  sphere  was  modeled  using  the  constitutive 
theory  of  Johnson  and  Cook  [9]  for  the  deviatoric 
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Fig.  7.  Ballistic  initial-boundary  value  problem:  (a)  coarse  mesh.  102  144  elements,  and  (b)  fine  mesh,  244  512  elements. 


Table  3 

Experimental  and  numerical  results,  ballistic  impact  and  fragmentation 


Experiment 

Simulation, 

coarse 

Simulation, 

fine 

Penetrator  Vr  (m/s) 

854 

814 

825 

Hole  diameter  (mm) 

21 

25 

26 

Crater  diameter  (mm) 

55 

52 

50 

Mass  loss  Mr  (kg) 

0.1 

0.134 

0.123 

Eroded  mass  M  (kg) 

— 

0.014 

0.014 

Avg.  frag,  velocity  v  (m/s) 

— 

133 

116 

mechanical  response,  with  a  Mie-Gruneisen  equation  of 
state  for  the  hydrostatic  response  [62].  As  the  tungsten 
sphere  did  not  fracture  or  degrade  significantly  during  the 
experiment,  failure  of  the  tungsten  was  suppressed  in  the 
simulations.  Properties  were  those  for  a  tungsten  alloy 
from  the  EPIC  library  [9,62],  apart  from  the  mass  density, 
which  was  set  to  exactly  match  that  of  the  experimental 
specimen.  It  is  noted  that  the  constitutive  updates  of  the 
material  responses  tended  to  require  little  execution  time 
relative  to  particle  search  algorithms  as  the  solutions 
proceeded  and  numerous  finite  elements  were  converted  to 
meshless  particles. 

Results  from  experiment  and  simulations  are  compared 
in  Table  3.  The  magnitude  of  the  residual  velocity  of  the 
spherical  projectile,  Fr,  agreed  within  5%  between 
simulations  and  experiment,  with  a  difference  of  less  than 
1.5%  between  simulations  with  fine  and  coarse  grids.  Hole 
and  crater  diameters  corresponded  to  the  exit  side  of  the 
target.  The  hole  diameter  was  estimated  experimentally  and 
in  the  simulations  as  the  visible  diameter  of  perforation. 
The  crater  diameter  was  measured  experimentally  as  is 
clear  from  visual  inspection  of  Fig.  8(a).  The  simulated 
crater  diameter  was  estimated  as  that  part  of  the  target 
surrounding  the  perforation  where  the  material  was  fully  or 
nearly  fully  damaged,  i.e.,  where  1.  This  region  is 
shown  by  darkened  grayscale  values  in  Fig.  8(b).  Con¬ 
verted  particles  are  deliberately  not  shown  in  Fig.  8(b).  It  is 


suggested  here  that  such  fully  damaged  material,  as  it 
supports  no  tensile  or  deviatoric  stresses,  would  simply  fall 
off  the  target  after  the  impact  event  due  to  the  force  of 
gravity,  leaving  behind  a  crater  in  the  concrete.  The  total 
mass  lost  in  the  experiment  was  determined  by  weighing 
the  target  before  and  after  the  test  and  subtracting  the  post¬ 
test  weight  from  the  initial  weight.  The  mass  loss  in  the 
simulation  was  computed  via 

Ml  =  J  pDdV,  (61) 

where  the  domain  of  integration  is  the  volume  of  target 
material,  and  damaged  material  was  considered  to 

contribute  to  mass  loss  following  the  above  arguments. 
The  eroded  mass  M  reported  in  the  simulations  was  simply 
the  summed  mass  of  all  particle  nodes,  while  the  average 
velocity  magnitude  of  the  fragment  cloud  was  computed  by 

0=  M~lS£2,mkvk,  (62) 

k 

where  the  velocity  if  of  each  particle  k  was  weighed  by  that 
particle’s  mass  nip.  Numerical  results  in  Table  3  correspond 
to  the  solution  time  t  =  150  ps,  as  values  of  these  quantities 
(e.g.,  average  projectile  and  fragment  velocities  and 

dimensions  derived  from  damage  contours)  did  not  change 
significantly  at  increments  beyond  this  instant  of  time  in 
the  calculations. 

Photographs  from  the  ballistic  experiment  are  shown  in 
Fig.  9.  The  projectile  is  exiting  the  target  and  traveling 
from  left  to  right  in  this  series  of  frames.  Initially,  some 
fragment  debris  is  propelled  ahead  of  the  projectile 
(Fig.  9(a)),  but  as  the  fragmentation  event  proceeds, 
the  projectile  passes  through  much  of  the  debris  cloud 
(Figs.  9(b)  and  (c)).  The  fragment  cloud  appears  close  to 
co-linear  with  the  path  of  the  projectile,  in  agreement 
with  prior  assumptions  in  the  statistical  model,  though 
there  is  some  spread  (i.e.,  non-zero  cone  angle).  The  fastest 
moving  particles  appear  small,  powder-like,  and  are 
thought  to  consist  mainly  of  pulverized  mortar.  Farger, 
slower-moving  particles  remain  closer  to  the  target  in  the 
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Fig.  8.  Post-event  crater:  (a)  experiment  [58]  and  (b)  simulation.  Darkened  areas  in  grid  (b)  denote  damage  state  variable  D. 


Fig.  9.  Fragmentation  event:  high-speed  photos  from  ballistic  experiment  [58], 


Fig.  9(c),  and  are  thought  to  consist  mainly  of  aggregate 
chunks.  Photographs  of  the  ballistic  event  were  not 
available  on  the  entrance  side  of  the  target,  though 
substantial  debris  was  observed  there  post-test. 

Results  from  a  simulation  are  visualized  in  Fig.  10,  at 
t  —  150  ps,  corresponding  to  the  experimental  image  in  Fig. 
9(c).  The  simulation  with  the  coarse  mesh  is  depicted; 
corresponding  results  for  the  fine  mesh  were  visually  very 
similar.  Concrete  particles  are  scaled  in  the  figure  by  the 
nominal  fragment  diameter  b,  as  computed  from  the 
energetic  theory,  Eq.  (32).  Particles  are  colored  by  velocity 
magnitude.  As  observed  in  the  experiment,  the  penetrator’s 
exit  velocity  exceeds  that  of  most  of  the  concrete  fragments. 
Larger  fragments  are  slower-moving  and  remain  close  to 
the  target,  while  fast  moving  fragments  are  relatively 
smaller  in  size.  Note  that  the  effects  of  aerodynamic 
resistance,  which  presumably  could  impart  drag  on  small 
particles  with  large  surface  area  to  volume  ratios,  are  not 
included  in  the  simulations.  Such  effects  would  presumably 
reduce  the  speed  of  very  fast  moving  particles  that  exceed 
the  projectile’s  velocity  in  the  simulation.  Contours  of 
damage  D  are  shown  in  the  concrete  target.  A  few  particles 
that  appear  stray  in  Fig.  10  were  in  fact  ejected  from  the 
sides  of  the  target  where  damage  propagated  to  the  edges. 
The  most  fully  damaged  material  logically  surrounds  the 
perforation,  though  the  damage  pattern  is  not  purely  axi- 
symmetric  due  to  the  non-symmetric  construction  of  the 
mesh  and  the  rectangular  geometry  of  the  target. 


Mass  probability  distributions  for  the  fragment  cloud 
are  shown  in  Fig.  11(a).  The  distributions  were  computed 
in  two  different  ways:  the  energetic  theory  of  Section  3.1 
and  the  statistical  physics-based  theory  of  Section  3.2.  For 
the  energetic  theory,  the  mass  of  all  fragments  comprising  a 
particle  k  was  computed  from  (32)  as 


Kq<x>c(D  -  T>t  )1  ' 
.  SKEpl/H2  \k 


(63) 


Subsequently,  distributions  were  generated  by  grouping 
fragments  into  bins  organized  by  mass: 


p(m\  <m<ni2 )  =  M  1  nij, 
i 


(64) 


where  j  spans  all  fragments  with  masses  in  the  range 
m\<m<m.2.  For  the  statistical  theory,  corresponding 
probability  distributions  were  found  from  (52)  as 


p(m\  < in Kmf) 


mp(m,  r)dr  dm 


“phr") 


m2 

mi 


(65) 


where  the  total  mass  M  is  listed  in  Table  3.  The  results  of 
Fig.  11(a)  may  be  interpreted  as  follows.  As  predicted  by 
the  energetic  theory  using  numerical  results  from  the 
coarse  grid,  for  example,  fragments  having  masses  between 
0.0001  and  0.001  kg  constitute  48%  of  the  total  mass  of  the 
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Fig.  10.  Fragmentation  event:  simulation  at  150ps.  Particle  diameter  scaled  by  fragment  dimension  and  particles  colored  by  velocity  magnitude. 
Penetrator  elements  colored  by  velocity  magnitude.  Concrete  target  elements  colored  by  damage  variable  D. 
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Fig.  11.  Predicted  fragment:  (a)  mass  and  (b)  velocity  distributions  at  150ps. 


fragment  cloud.  Slight  differences  in  mass  probabilities 
among  simulations  with  coarse  and  fine  grids  are  evident  in 
the  energetic  theory,  on  the  order  of  a  few  percent  at  most 
for  any  particular  bin.  However,  since  an  identical  value  of 
M  was  computed  in  both  simulations,  the  statistical  theory 
yielded  identical  results  for  coarse  and  fine  discretizations. 
Note  that  for  all  four  cases  (i.e.,  two  theories  and  two  grid 
sizes),  the  bin  containing  the  largest  mass  fraction  of 
fragments  was  that  for  which  fragments  spanned  the  range 
of  0.0001  kg  <m<  0.001  kg. 


Velocity  probability  distributions  of  the  fragment  debris 
are  shown  in  Fig.  1 1(b).  These  are  mass-weighted  distribu¬ 
tions.  For  the  energetic  theory,  the  probability  was 
computed  by 

p(v i  <v< V2)  =  M~l  Y mk,  (66) 

k 

where  the  index  k  here  identifies  a  particle  whose  velocity  is 
within  the  range  v\<v<V2,  with  mass  of  that  particle 
denoted  by  m^.  For  the  statistical  physics  theory,  individual 
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particle  velocities  were  not  considered  directly.  Instead  the 
global  center  of  mass  velocity  v  of  Eq.  (62)  and  Table  3  was 
substituted  into  (57)  to  compute  the  kinetic  energy 
distributed  among  the  fragments.  The  probability  was  then 
estimated  from  Eq.  (52)  as 

/*  v>2  r  oo 

p(v\<v<V2)  —  N~l  /  /  p(m,  v)dmdv  —  ^(3N)/(4E) 

J  v  i  JO 

X  [v/  f  {N / M)  +  (3Nv2)/(4E)]  | £ .  (67) 

Fig.  1 1(b)  shows  the  predicted  mass  fraction  of  fragments 
within  a  particular  velocity  range.  For  example,  from  the 
energetic  theory  with  a  coarse  mesh,  47%  of  the  mass  of  the 
cloud  should  exhibit  a  velocity  magnitude  in  the  range 
0-100  m/s,  and  34%  in  the  range  of  100-200  m/s.  Slight 
differences  between  velocity  v  for  simulations  with  coarse 
and  fine  meshes  led  to  differences  in  their  statistical  physics- 
based  predictions,  for  example  on  the  order  of  a  4%  for  the 
computed  value  of  p  for  the  mass  fraction  of  fragments  with 
r<100m/s.  Predictions  from  the  energetic  theory  differed 
slightly  more  with  mesh  density,  on  order  of  8%,  for 
example,  for  the  mass  fraction  of  fragments  with 
r<100m/s.  Note  that  for  all  four  cases,  the  bin  of  range 
0<y<100m/s  was  predicted  to  contain  the  largest  mass 
fraction  of  fragments. 

5.2.  Dynamic  impact  crushing 

Because  quantitative  data  on  fragment  mass  distribu¬ 
tions  were  not  available  from  the  experiment  [58]  simulated 
above,  a  second  simulation  was  performed  to  further  judge 
the  accuracy  and  limitations  of  aspects  of  the  present 
model.  The  experiment  was  conducted  elsewhere,  as 
described  in  [3],  in  an  investigation  of  methods  for  reducing 
concrete  rubble  to  recyclable  form  via  impact  crushing.  In 
the  experiment  [3],  a  concrete  sphere  was  propelled  at 
moderate  velocity  into  a  target  plate  within  a  cylindrical 
catch  tank;  upon  impact  of  the  sphere  with  the  plate,  the 
concrete  fractured  and  fragmented  into  debris  that  were 
collected  from  the  tank.  Multiple  passes  through  the 
apparatus  were  used  to  further  reduce  the  rubble.  Here, 
only  the  first  pass  is  simulated,  and  the  interactions 
between  concrete  projectile  and  metallic  tank  are  limited 
to  the  initial  impact  of  the  sphere  with  a  flat  plate. 
Subsequent  interactions  that  would  presumably  occur  in 
the  experiment,  such  as  impact  of  the  fragments  with  the 
sides  of  the  tank  upon  ricochet,  are  not  modeled  here. 

The  same  modeling  procedures  developed  in  Sections 
2-4  were  used  here.  Elastic,  plastic,  and  damage  evolution 
relationships  and  properties  for  SAC-7  concrete  [6]  were 
employed  in  the  simulation,  as  experimental  data  on  the 
strength  of  the  concrete  tested  in  [3]  were  not  reported.  For 
this  reason,  discrepancies  may  naturally  emerge  between 
results  of  simulation  and  experiment.  Specifically,  it  is 
noted  that  the  compositions  appear  different  in  [3]  and  [6], 
with  different  types  of  aggregate  and  mortar.  In  order  to 
predict  the  fragment  mass  distributions,  only  the  method 


predicated  upon  a  local  energy  balance  (Section  3.1)  was 
used  here.  Since  the  fragment  clouds  are  far  from  ID  (52) 
or  spherical  [57],  analytical  statistical  physics-based  the¬ 
ories  of  the  sort  derived  in  Section  3.2  do  not  immediately 
apply. 

The  problem  consisted  of  a  spherical  concrete  projectile 
of  diameter  150nnn  impacted  against  a  steel  wall  of 
dimensions  100  mm  x  750  mm  x  750  mm  at  a  normal 
velocity  of  55  m/s.  Two  meshes  of  increasingly  fine 
discretizaton  were  used.  The  coarse  mesh  consisted  of 
34  272  elements,  while  the  fine  mesh  included  8 1 408 
elements.  As  in  the  previous  problem  of  Section  5.1,  all 
elements  were  composite  tetrahedral  [53],  and  mesh 
densities  in  the  target  were  reduced  gradually  with 
increasing  distance  from  the  location  of  impact.  Friction¬ 
less  contact  between  projectile  and  target  was  assumed,  and 
adiabatic  conditions  were  invoked,  with  a  uniform  initial 
temperature  of  294  K.  The  target  plate  was  assumed  to 
consist  of  4340  steel,  addressed  here  with  standard 
Johnson-Cook  plasticity  [9]  and  fracture  [10]  models  and 
properties  from  the  EPIC  material  library.  The  coarse 
mesh  is  shown  in  Fig.  12(a). 

The  deformed  coarse  mesh,  with  converted  particles,  is 
shown  in  Fig.  12(b),  at  t  =  5  ms  after  impact.  Fragment 
mass  distributions  remained  nearly  static  at  later  solution 
times.  Particle  diameters  are  scaled  by  their  corresponding 
fragment  size  b  in  Fig.  12(b).  Note  that  each  visually 
smaller  particle  in  fact  often  represents  a  cloud  of  many 
small  particles,  since  the  nodal  mass  of  each  particle  is 
distributed  among  the  multiple  fragments  corresponding  to 
that  particle.  From  Fig.  12(b)  a  mixed  distribution  of 
fragment  sizes  is  apparent.  The  concrete  debris  consists  of 
small  faster-moving  fragments  on  the  periphery,  large 
generally  slower-moving  fragments  near  the  center  of  the 
target  plate,  and  several  intact  pieces  of  the  projectile.  The 
latter  are  finite  elements  that  did  not  fail  or  convert  to 
particles  in  the  simulation. 

Predicted  fragment  mass  distributions,  computed  from 
Eq.  (64),  are  compared  with  those  obtained  from  the 
experiment  [3]  in  Fig.  13.  Bins  used  to  group  the  fragments 
here  are  the  same  used  in  [3],  where  fragment  diameters 
from  [3]  have  been  converted  to  masses  in  Fig.  13  via  the 
mass  density  of  SAC-7  concrete  listed  in  Table  1, 
p0  =  2440  kg/m3.  From  Fig.  13,  the  model  appears  to 
predict  the  fragment  mass  distributions  with  reasonable 
accuracy.  The  relatively  large  concentration  p  ss  0.4  of 
larger  fragments  in  the  range  of  0.02  kg</«<  1.3  kg  results 
mainly  from  the  unconverted  finite  elements  of  the  concrete 
sphere  (perhaps  damaged  but  not  failed).  These  pieces  were 
included  in  the  fragment  count,  since  in  the  experiment  the 
entire  specimen  was  apparently  sieved  during  fragment 
collection  from  the  catch  tank.  The  relatively  large 
concentration  p  ^  0.35  in  the  intermediate  range  of 
0.00004  kg <m< 0.0025  kg  emerges  from  the  balance  of 
fracture  energy  and  kinetic  energy  fueling  the  local 
fragmentation  process  in  (32);  i.e.,  most  of  the  converted 
particles  fall  into  this  bin.  The  fraction  of  smaller 
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Fig.  12.  Impact  crushing  of  concrete  sphere:  (a)  initial  conditions  and  (b)  fragmentation  event  at  5  ms,  with  particle  diameters  scaled  by  fragment 
dimension. 
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Fig.  13.  Fragment  mass  distributions  from  impact  crushing  of  concrete 
sphere  (experimental  data  from  [3]). 

fragments,  «7<5.2(10)~6  kg,  is  vastly  under-predicted  by 
the  theory.  This  is  because  local  strain  rates  achieved 
during  integration  cycles  over  which  D>Dj  were  of 
insufficient  magnitude  to  generate  very  small  fragments 
via  (32).  This  effect  highlights  the  limitation  of  the  present 
modeling  approach:  only  dynamic  high-rate  fragmentation 
problems  can  be  simulated,  as  the  local  kinetic  energy  of 
the  material,  balanced  by  the  energy  dissipated  by  micro¬ 
cracking,  provides  the  sole  impetus  for  break  up  into 
numerous  fragments.  It  is  noted  that  very  small  fragments 
could  feasibly  be  generated  in  the  experiments  by  interac¬ 
tions  with  the  walls  of  the  catch  tank  upon  ricochet  and 
settling,  effects  not  modeled  here.  It  is  also  noted  that 
explicit  radial  crack  patterns  reported  from  high-speed 
photographs  in  [3]  cannot  be  reproduced  by  the  present 
framework  based  on  continuum  damage  mechanics.  Some 
mesh  dependence  is  evident,  but  not  excessive,  in  Fig.  13: 
differences  among  predictions  of  p(m)  obtained  from  fine 


and  coarse  meshes  of  no  more  than  6%  are  apparent  in 
any  bin. 

Fragment  velocity  data  were  not  available  from  the 
experimental  description  [3];  hence  such  data,  though 
available  from  the  computation,  are  not  presented  here. 
Experimental  velocity  data  appear  scarce  for  dynamic 
fragmentation  of  pure  (not  reinforced)  concrete  of  interest 
here.  Flowever,  some  experimental  data  are  available  for 
colliding  rocks  [57,71]  and  bar-reinforced  concrete  walls 
undergoing  explosive  loading  [2]. 

6.  Conclusions 

A  constitutive  theory  for  deforming  and  fragmenting 
crushable  solids  has  been  developed.  The  theory  includes 
the  following  novel  combination  of  features:  finite  defor¬ 
mation  kinematics  of  elasticity,  plasticity,  and  inelastic 
compression;  an  energy  density  function  with  contributions 
from  nonlinear  elasticity,  damage,  and  porosity;  and 
kinetic  relations  for  evolution  of  inelastic  deformation, 
pore  compression,  and  damage  satisfying  both  the 
energy  balance  and  entropy  inequality.  The  theory  also 
features  two  alternative,  physically  motivated  methods  to 
compute  fragment  size  and  velocity  distributions.  In  the 
energetic  method,  a  local  balance  between  expansion 
kinetic  energy  and  fracture  energy  associated  with  damage 
evolution  leads  to  a  local  fragment  size.  In  the  statistical 
method,  global  entropy  maximization  leads  to  a  joint 
mass-velocity  distribution  function.  The  theory  has  been 
implemented  in  a  Lagrangian  FE  setting  with  GPA, 
whereby  highly  distorted  elements  are  converted  to 
interacting  particle  nodes.  In  the  implementation  of  the 
plasticity  theory,  small  elastic  strain  assumptions  are  used, 
though  the  hydrostatic  response,  computed  separately,  is 
considered  valid  to  large  pressures.  Local/global  solution 
data  from  converted  GPA  nodes  are  used  to  compute 
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fragment  size  and  velocity  distributions  for  the  energetic/ 
statistical  methods,  respectively. 

The  model  was  first  used  to  study  the  impact  of  a  small 
tungsten  sphere  into  a  thin  concrete  target.  Residual 
projectile  velocities,  hole  and  crater  dimensions,  and  target 
mass  loss  were  comparable  between  experiment  and  two 
simulations  of  different  mesh  density.  Attributes  of  the 
fragment  debris  clouds  agreed  qualitatively  among  simula¬ 
tions  and  experiment,  with  the  fastest  particles  usually 
being  smallest  in  size.  Distributions  of  mass  and  velocity 
were  computed  using  the  two  thermodynamically  moti¬ 
vated  methods  (energetic  and  statistical),  and  two  grid 
spacings  (coarse  and  fine).  Mass-weighted  distributions 
were  qualitatively  similar  among  all  cases,  with  the  bin 
corresponding  to  the  largest  fraction  of  fragments  having  a 
particular  mass  or  velocity  range  the  same  regardless  of 
method  or  mesh  size.  Predictions  of  mass  distribution  were 
fairly  mesh-independent,  with  velocity  slightly  more  so 
(e.g.,  maximum  velocity  probability  differences  on  the 
order  of  8%  for  a  140%  increase  in  number  of  elements).  It 
remains  to  be  seen  which  fragmentation  theory  most 
closely  replicates  live  ballistic  perforation  experiments  as 
experimental  mass  and  velocity  distribution  data  become 
more  available. 

In  a  second  simulation,  the  energetic  method  was  used  to 
predict  fragment  mass  distributions  upon  impact  crushing 
of  a  concrete  sphere.  Following  quantitative  comparison 
with  experimental  data,  the  theory  predicted  mass  dis¬ 
tributions  of  larger  fragments  with  reasonable  success, 
despite  uncertainty  in  the  strength  properties  of  the 
concrete  tested  in  the  experiment.  Flowever,  the  model 
failed  to  capture  very  small  fragment  debris  produced  in 
the  experiment,  as  sufficiently  high  deformation  rates 
needed  to  yield  very  small  fragments  were  not  attained  in 
the  simulation.  This  artifact  of  the  model  illustrates  that 
modifications  to  the  present  framework  are  needed  to 
properly  address  failure  and  fragmentation  under  quasi¬ 
static  and  low-strain-rate  loading  conditions. 
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